2

I have a NetCDF file that contains a series of HDF5 bands. I am looking to project one of these bands contained within this file. I've attempted to use the following command to achieve this:

gdalwarp -geoloc -of GTIFF -t_srs EPSG:4269 HDF5:"zz_ssmv01020S_latest.nc"://Surface_temperature

I get the following error: 

ERROR 1: Unable to compute a GEOLOC_ARRAY based transformation between pixel/line and georeferenced coordinates for HDF5:zz_ssmv01020S_latest.nc://Surface_temperature

My goal is to project this grid into the North American Datum reference (EPSG:4269), as metadata already indicates the coordinates are in lat lon. Any ideas as to why gdalwarp fails to project this?

The readout from the desired band to extract using gdalinfo is as follows: 

    Size is 8192, 4096
    Coordinate System is `'
    Metadata:
      Conventions=CF-1.6
      created_date=2019-03-26 16:02:56
      institution=National Operational Hydrologic Remote Sensing Center, Chanhassen, MN
      references=http://www.nohrsc.noaa.gov/technology/
      source=National Snow Analyses
      _nc3_strict=1
    Corner Coordinates:
    Upper Left  (    0.0,    0.0)
    Lower Left  (    0.0, 4096.0)
    Upper Right ( 8192.0,    0.0)
    Lower Right ( 8192.0, 4096.0)
    Center      ( 4096.0, 2048.0)
    Band 1 Block=8192x1 Type=Int16, ColorInterp=Undefined
      Metadata:
        Surface_temperature_coordinates=lat lon
        Surface_temperature_grid_mapping=crs
        Surface_temperature_long_name=Surface temperature
        Surface_temperature_maximum_data_value=299.6
        Surface_temperature_minimum_data_value=238.9
        Surface_temperature_missing_value=-32768
        Surface_temperature_scale_factor=0.01
        Surface_temperature_standard_name=air_temperature
        Surface_temperature_start_date=2019-03-26 11:00:00
        Surface_temperature_stop_date=2019-03-26 11:00:00
        Surface_temperature_title=Scaled corrected Surface temperature analysis
        Surface_temperature_units=Kelvins
        Surface_temperature__FillValue=-32768

The gdalinfo result from the entire file is as follows: 

Driver: HDF5/Hierarchical Data Format Release 5
Files: zz_ssmv01020S_latest.nc
Size is 512, 512
Coordinate System is `'
Metadata:
  Conventions=CF-1.6
  created_date=2019-03-26 16:02:56
  crs_grid_mapping_name=latitude_longitude
  crs_horizontal_datum=WGS84
  crs_inverse_flattening=298.257222096042
  crs_longitude_of_prime_meridian=0
  crs_semi_major_axis=6378137
  crs_semi_minor_axis=6356752.31414
  institution=National Operational Hydrologic Remote Sensing Center, Chanhassen, MN
  lat_axis=Y
  lat_bounds=lat_bounds
  lat_CLASS=DIMENSION_SCALE
  lat_comment=location of pixel center
  lat_long_name=latitude
  lat_NAME=lat
  lat_origin_offset=0.00416666666667146
  lat_REFERENCE_LIST=
  lat_resolution=0.00833333333333333
  lat_standard_name=latitude
  lat_units=degrees_north
  lat__Netcdf4Dimid=0
  lon_axis=X
  lon_bounds=lon_bounds
  lon_CLASS=DIMENSION_SCALE
  lon_comment=location of pixel center
  lon_long_name=longitude
  lon_NAME=lon
  lon_origin_offset=0.00416666666666052
  lon_REFERENCE_LIST=
  lon_resolution=0.00833333333333333
  lon_standard_name=longitude
  lon_units=degrees_east
  lon__Netcdf4Dimid=1
  nv_CLASS=DIMENSION_SCALE
  nv_NAME=This is a netCDF dimension but not a netCDF variable.         2
  nv_REFERENCE_LIST=
  nv__Netcdf4Dimid=3
  references=http://www.nohrsc.noaa.gov/technology/
  source=National Snow Analyses
  Surface_temperature_coordinates=lat lon
  Surface_temperature_grid_mapping=crs
  Surface_temperature_long_name=Surface temperature
  Surface_temperature_maximum_data_value=299.6
  Surface_temperature_minimum_data_value=238.9
  Surface_temperature_missing_value=-32768
  Surface_temperature_scale_factor=0.01
  Surface_temperature_standard_name=air_temperature
  Surface_temperature_start_date=2019-03-26 11:00:00
  Surface_temperature_stop_date=2019-03-26 11:00:00
  Surface_temperature_title=Scaled corrected Surface temperature analysis
  Surface_temperature_units=Kelvins
  Surface_temperature__FillValue=-32768
  time_CLASS=DIMENSION_SCALE
  time_long_name=time
  time_NAME=time
  time_REFERENCE_LIST=
  time_units=hours since 2019-03-26 11:00:00
  time__Netcdf4Dimid=2
  _nc3_strict=1
Subdatasets:
  SUBDATASET_1_NAME=HDF5:"zz_ssmv01020S_latest.nc"://Surface_temperature
  SUBDATASET_1_DESC=[1x4096x8192] //Surface_temperature (16-bit integer)
  SUBDATASET_2_NAME=HDF5:"zz_ssmv01020S_latest.nc"://lat_bounds
  SUBDATASET_2_DESC=[4096x2] //lat_bounds (64-bit floating-point)
  SUBDATASET_3_NAME=HDF5:"zz_ssmv01020S_latest.nc"://lon_bounds
  SUBDATASET_3_DESC=[8192x2] //lon_bounds (64-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)
Improve this question
2

1 Answer 1

Reset to default
1

Unfortunately, gdalwarp -geoloc still does not recognize lat and lon variables by itself (or at least I have not yet found out how to tell it ...).

But you can work around that using an intermediate VRT with <Metadata domain="GEOLOCATION"> field.

  1. Export the dataset as not georeferenced GTiff:

    nodata_value=-32768
scale_factor=0.01

gdal_calc.py --format=GTiff --type=Float32 \
    -A "HDF5:zz_ssmv01020S_latest.nc://Surface_temperature" \
    --calc="(A!=$nodata_value)*A*$scale_factor+(A==$nodata_value)*A" \
    --outfile tmpfile1.tif

  2. Convert the result to VRT

    gdal_translate -of VRT tmpfile1.tif tmpfile2.vrt

  3. Add GEOLOCATION information to VRT:

    <Metadata domain="GEOLOCATION">
  <MDI key="X_DATASET">HDF5:zz_ssmv01020S_latest.nc://lon</MDI>
  <MDI key="X_BAND">1</MDI>
  <MDI key="Y_DATASET">HDF5:zz_ssmv01020S_latest.nc://lat</MDI>
  <MDI key="Y_BAND">1</MDI>
  <MDI key="PIXEL_OFFSET">0</MDI>
  <MDI key="LINE_OFFSET">0</MDI>
  <MDI key="PIXEL_STEP">1</MDI>
  <MDI key="LINE_STEP">1</MDI>
</Metadata>

    (I do this automated using sed -i '2i...' tmpfile2.vrt)

  4. Run gdalwarp -geoloc on the VRT

    gdalwarp -of GTiff -geoloc -s_srs EPSG:4326 -t_srs EPSG:4269 \
  -dstnodata "$nodata_value" \
  tmpfile2.vrt result.tif
Improve this answer
1
  • Thanks for writing out a detailed process @RoVo . The code works well, however at step 4, I come to the following error: ERROR 1: X_BAND and Y_BAND do not have the same dimensions. Any insight as to what's throwing this error?
    – TornadoEric
    Commented Mar 27, 2019 at 14:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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