Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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 `'
Metadata:
  NC_GLOBAL#Conventions=CF-1.4
  NC_GLOBAL#Satellite Sensor=G-13 IMG    
  NC_GLOBAL#Source=McIDAS Area File
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"goes13.2013.100.165517.BAND_04.nc":auditTrail
  SUBDATASET_1_DESC=[3x80] auditTrail (8-bit character)
  SUBDATASET_2_NAME=NETCDF:"goes13.2013.100.165517.BAND_04.nc":data
  SUBDATASET_2_DESC=[1x665x2036] data (32-bit floating-point)
  SUBDATASET_3_NAME=NETCDF:"goes13.2013.100.165517.BAND_04.nc":lat
  SUBDATASET_3_DESC=[665x2036] lat (32-bit floating-point)
  SUBDATASET_4_NAME=NETCDF:"goes13.2013.100.165517.BAND_04.nc":lon
  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 `'
Metadata:
  data#coordinates=lon lat
  data#long_name=0-255 Brightness Temperature
  data#type=VISR
  NC_GLOBAL#Conventions=CF-1.4
  NC_GLOBAL#Satellite Sensor=G-13 IMG    
  NC_GLOBAL#Source=McIDAS Area File
  NETCDF_DIM_EXTRA={time}
  NETCDF_DIM_time_DEF={1,4}
  NETCDF_DIM_time_VALUES=1365615900
  time#long_name=seconds since 1970-1-1 0:0:0
  time#units=seconds since 1970-1-1 0:0:0
Geolocation:
  LINE_OFFSET=0
  LINE_STEP=1
  PIXEL_OFFSET=0
  PIXEL_STEP=1
  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"]]
  X_BAND=1
  X_DATASET=NETCDF:"goes13.2013.100.174518.BAND_04.nc":lon
  Y_BAND=1
  Y_DATASET=NETCDF:"goes13.2013.100.174518.BAND_04.nc":lat
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
  Metadata:
    coordinates=lon lat
    long_name=0-255 Brightness Temperature
    NETCDF_DIM_time=1365615900
    NETCDF_VARNAME=data
    type=VISR

Any help will be greatly appreciated!

share|improve this question
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? –  kyle Jul 23 '13 at 14:56

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.