I am trying to remap from Geostationary to Lambert conformal using gdalwarp. My input data is in netcdf, and is in geographic coordinates (degrees) and I would like to output the remapped data to netcdf. I have created a corresponding vrt file for the input netcdf data. Gdalwarp will output the netcdf file, but the output data is all zeros and I receive the following error:

Creating output file that is 5120P x 5120L.
Processing input file netcdf.vrt.
ERROR 1: Too many points (441 out of 441) failed to transform,
unable to compute output bounds.
Warning 1: Unable to compute source region for output window 0,0,5120,5120, skipping.
0...10...20...30...40...50...60...70...80...90...100 - done.

I attempted the following command:

/usr/bin/gdalwarp -s_srs "+proj=geos +h=35785831 +lon_0=-75 +x_0=-0.151844 +y_0=0.151844 +a=6378140 +b=6356754.99999591 +units=degrees +no_defs" -t_srs "+proj=lcc +ellps=clrk66 +a=6378137 +b=6378137 +e=0.0818191910435 +lat_0=24.9999 +lon_0=-95 +lat_1=24.9999 +lat_ts=25.0001 +units=meters +no_defs" -te -1952976.3246 -828316.5944 3248431.6754 4373091.4056 -of netCDF -geoloc -overwrite -r bilinear -ts 5120 5120 netcdf.vrt out.nc

Can gdalwarp remap from geographic coordinates to projected? Or do I need to translate geographic to projected first? Also, can gdalwarp read projection information straight from netcdf or do you NEED to write to .vrt first?

Here is what gdalinfo outputs from the input file: (it is a GOES 13 file from CLASS)

Size is 512, 512
Coordinate System is `'
  NC_GLOBAL#Satellite Sensor=G-13 IMG    
  NC_GLOBAL#Source=McIDAS Area File
  SUBDATASET_1_DESC=[3x80] auditTrail (8-bit character)
  SUBDATASET_2_DESC=[1x665x2036] data (32-bit floating-point)
  SUBDATASET_3_DESC=[665x2036] lat (32-bit floating-point)
  SUBDATASET_4_DESC=[665x2036] lon (32-bit floating-point)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

And additional gdal info on the data variable:

Driver: netCDF/Network Common Data Format
Files: goes13.2013.100.174518.BAND_04.nc
Size is 2036, 665
Coordinate System is `'
  data#coordinates=lon lat
  data#long_name=0-255 Brightness Temperature
  NC_GLOBAL#Satellite Sensor=G-13 IMG    
  NC_GLOBAL#Source=McIDAS Area File
  time#long_name=seconds since 1970-1-1 0:0:0
  time#units=seconds since 1970-1-1 0:0:0
  SRS=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  665.0)
Upper Right ( 2036.0,    0.0)
Lower Right ( 2036.0,  665.0)
Center      ( 1018.0,  332.5)
Band 1 Block=2036x1 Type=Float32, ColorInterp=Undefined
  NoData Value=9.96920996838686905e+36
    coordinates=lon lat
    long_name=0-255 Brightness Temperature

Any help will be greatly appreciated!

  • 1
    The geos projection won't use degrees; try meters. Where are you getting the +x_0/+y_0 values? Based on gdalinfo, I'm not sure the input raster is georeferenced at all. In the target srs, you have +a = +b which is a sphere, but also set +e. However, +ellps is for an entirely different ellipsoid. The various +lat values all seem odd too. lat_ts is lat of true scale so a standard parallel just like lat_1. – mkennedy Jul 22 '13 at 15:54
  • Thank you. I will try using meters. I am getting the x_0 and y_0 (scale and offsets) from the GOES definition although these are not required inputs for +proj=geos, so I can try taking them out. And thanks for pointing out the +e ellipsoid addition The lat definitions for the t_srs are for AWIPS definiton of lambert (a specific output size). I'll add what gdal info spits out for the specific data variable to the question post (too long for comment) – Katie J Jul 22 '13 at 16:08
  • The AWIPS definition I'm referring to is described on this page: nws.noaa.gov/noaaport/html/icdtb48_2.html (the first one is the Lambert I am attempting to remap to) – Katie J Jul 22 '13 at 16:35
  • 1
    Hmmm, so it says lat/lon WGS84, but the reported corner coordinates concern me because they're just raw cell values. The LCC is a tangent case--single standard parallel/latitude of origin are all at 25N. I haven't worked with any of this data, so I'm just going by the metadata info. – mkennedy Jul 22 '13 at 18:06
  • The image is not georeferenced, but a source srs is supplies. A few questions: * Can you run with CPL_DEBUG=GDAL_netCDF? So CPL_DEBUG=GDAL_netCDF /usr/bin/gdalwarp ... I suspect that there may be a problem with the geolocation arrays. * Can you make your data available? – user10353 Jul 23 '13 at 14:56

If your datasource has longitude and latitude cell values as separate subdatasets, try building a vrt file manually to reproject, as described in

Unable to warp HDF5 files and gdal for reprojection of VSCMO VIIRS data

In a first step, reproject to EPSG:4326 using the subdatasets, then to whatever CRS you want.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.