Ed,
Glad to hear you're making good progress. I wasn't clear in describing
the proj call. In the GeoTransform, (X0,Y0) is the upper left corner of
the source grid in coordinates of the source projection. From the
lon/lat grids from your source dataset, you can find the lon/lat of the
upper left corner, (lon0,lat0), and lower right corner (lon1,lat1) of
your grid. The question is "what coordinates, (X0,Y0) and (X1,Y1), in
the source projection correspond to (lon0,lat0) and (lon1,lat1)?" The
proj call below answers questions of this form:
proj +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000
+b=6370000 <<EOF
LON0 LAT0
LON1 LAT1
EOF
Where you should replace LON0, LAT0, LON1, LAT1 with the
upper_left_longitude, upper_left_latitude, lower_right_longitude, and
lower_right_latitude from your grid, respectively.
Good luck.
Steve
Ed Fialkowski wrote:
Steve, that write up was excellent. Thank you for that. You
definitely weren't telling me anything I already know.
If I could just get a little more help, I think I'm so close. My
output file seems to have problems. The resulting image is 25k, but
it's just black.
Gdalinfo on the resulting geotiff (test.tif):
Driver: GTiff/GeoTIFF
Size is 119, 152
Coordinate System is:
GEOGCS["WGS 84",
DATUM["unknown",
SPHEROID["unnamed",6378137,298.2572235630016]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]]
Origin = (-113.117651610910784,60.203870178822328)
Pixel Size = (0.135957945826561,-0.135957945826561)
Metadata:
AREA_OR_POINT=Area
Corner Coordinates:
Upper Left (-113.1176516, 60.2038702) (113d 7'3.55"W, 60d12'13.93"N)
Lower Left (-113.1176516, 39.5382624) (113d 7'3.55"W, 39d32'17.74"N)
Upper Right ( -96.9386561, 60.2038702) ( 96d56'19.16"W, 60d12'13.93"N)
Lower Right ( -96.9386561, 39.5382624) ( 96d56'19.16"W, 39d32'17.74"N)
Center (-105.0281538, 49.8710663) (105d 1'41.35"W, 49d52'15.84"N)
Band 1 Block=119x68 Type=Byte, ColorInterp=Gray
(Coordinates look pretty good, just about over US)
I wasn't completely clear on the correct syntax for the proj call, the
first few I tried gave me a 0 dx and dy, so here's the best I got:
proj +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000
+b=6370000 <<EOF
0_0_0 0_0_0
148_ny_nx 112_ny_nx
EOF
replaced ("LON" with 148 and "LAT" with 112, my lower right point is
148.0,112.0)
which gives me:
10908379.76 1690481.54_0_0 0_0_0
-11877294.73 3954472.04_ny_nx 112_ny_nx
gives me dx=-6591.2583, dy=20396.3108
You said geoTransform is x0-dx/2, dx, 0, y0-dy/2, dy. I'm hoping x0
and y0 are 0 (upper left?)
<GeoTransform>3295.62915, -6591.2583, 0.0, -10198.1554, 0.0,
20396.3108</GeoTransform>
in form of x0-dx/2, dx, 0, y0-dy/2, 0, dy
I hate to ask you to check all my syntax after such amazing help, but
I would be incredibly grateful.
Thanks again.
-Ed
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
