[Gdal-dev] problems warping a netcdf file

[Steve Gaffigan]

Ed,

Looking at your gdalwarp call, it seems that you've omitted the 
projection parameters.  Based on your description of the projection, I 
would expect a gdalwarp call such as the following:

gdalwarp -s_srs "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 
+a=6370000 +b=6370000 +x_0=0 +y_0=0" -t_srs "+proj=longlat" 
NETCDF:"test_0910.nc":PM25 test_out.nc

The real problem, however, is related to the GeoTransform of your source 
dataset.  Though you didn't show the full gdalinfo output, based on your 
description of the dataset I expect gdalinfo shows the following under 
"Corner Coordinates":

Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  112.0)
Upper Right (  148.0,    0.0)
Lower Right (  148.0,  112.0)
Center      (  74.0,  56.0)

If this is the case, then it indicates that you'll need to additionally 
provide gdalwarp with the GeoTransform of your dataset, as the above 
shows that gdal doesn't know anything more about your grid than 
column/row (e.g. pixel/line) index values.  There are a number of ways 
to specify both projection and geotransform in the netcdf file itself, 
but maybe the easiest thing at this point would be to create, in the 
directory containing your file, a VRT definition file containing the 
following:

  <VRTDataset rasterXSize="148" rasterYSize="112">
   <SRS>+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 
+b=6370000</SRS>
   <GeoTransform>X0, DX, 0.0, Y0, 0.0, DY</GeoTransform>
   <VRTRasterBand>
     <SimpleSource>
       <SourceFilename 
relativeToVRT="1">NETCDF:test_0910.nc:PM25</SourceFilename>
     </SimpleSource>
   </VRTRasterBand>
</VRTDataset>

In the GeoTransform tag, you'll need to replace X0,DX,Y0,DY with the 
correct values for your dataset.  Sorry if I'm telling you what you 
already know, but X0 and Y0 are the upper left (line=0,column=0) easting 
and northing coordinates of your grid, respectively.  dx and dy are the 
spacings between columns and lines, respectively.  All these parameters 
are in units of meters relative to the native projection of your 
dataset.  You might find these values in file or source metadata. 
Otherwise, if you have the latitude and longitude grids for your source 
dataset, you can calculate the required parameters using known 
projection transformation equations.  There are a number of ways to do 
this, but if you have the proj4 library on your system you can do it on 
the command line as follows:

proj +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 
+b=6370000 <<EOF
LON_0_0 LAT_0_0
LON_ny_nx LAT_ny_nx
EOF

Where you'll need to replace these LON,LAT values with the relevant 
upper left (0,0) and lower right (ny,nx) values.  The output will be two 
lines showing easting,northing values for the corners:

X0 Y0
X1 Y1

 From this you can calculate your dx and dy:

dx=(X1-X0)/147
dy=(Y1-Y0)/111

Then the GeoTransform is:
x0-dx/2, dx, 0, y0-dy/2, dy

Finally, your new gdalwarp call would be:

gdalwarp -t_srs "+proj=longlat" test_0910.vrt test_out.nc

Steve


edfialk wrote:
> Hello, I'm fairly new to reprojection and my experience with gdal is limited,
> so I'll do my best to keep up.  I have tried searching through a number of
> messages without finding anything quite like my problem.
> 
> I've been completely unsuccessful trying to reproject a netcdf file for some
> time now and I'm becoming fairly desperate for some help.  I'm going for
> changing "Lambert Conformal with a perfect sphere geoid with earth
> radius=6370km centered at 40N, 97W with standard latitudes at 33N, 45N to
> geographic lat, long (flat) projection."
> 
> First, a little gdal info on one particular subdataset for the netcdf file:
> Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
> Driver: netCDF/Network Common Data Format
> Size is 148, 112
> Coordinate System is `'
> ...
> ad then coordinates, cdate, ctime, xcent, ycent, etc. etc.  I was wondering
> if NC_Global or the coordinate system being `' might cause a problem.  
> 
> and then, running the command (going to mercator just to get a successful
> reprojection, I'll deal with lat/lon later):
> gdalwarp -s_srs "+proj=lcc" -t_srs "+proj=merc" NETCDF:"test_0910.nc":PM25
> NETCDF:"test_out.nc":PM25
> 
> gives me:
> Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
> ERROR 1: Unable to compute a transformation between pixel/line
> and georeferenced coordinates for NETCDF:test_0910.nc:PM25.
> There is no affine transformation and no GCPs.
> 
> I'm fairly stuck.  If anyone has any ideas at all, I'd really appreciate it.
> Thanks!
> -Ed