[Gdal-dev] problems warping a netcdf file

[Steve Gaffigan]

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