1

I'm using the python API to GDAL to try to do a simple geographic warp of a netCDF4 satellite file to a small sub region. My python program fails to produce any output but does not throw any exceptions or errors. Any ideas?

The command line operation of gdalwarp works fine and produces meaningful netCDF output (for a single band or variable). But the python version does not. THe output 'ds' is "None' and no netCDF output file is created. Here is a snapshot:

    from osgeo import gdal, gdalconst


kwargs = {'format': 'netCDF', 'copyMetadata': True, 
          'outputBounds': '[-29 -2 -25 2]', 
          'dstSRS':'+proj=longlat +datum=WGS84 +no_defs'
         }
          

src_filename = './subset_20180729030000-JPL-L2P_GHRSST-SSTskin-MODIS_A-N-v02.0-fv01.0.nc'
nc_file = 'NETCDF:' + src_filename
src = gdal.Open(nc_file, gdalconst.GA_ReadOnly)

subDatasets = src.GetSubDatasets()
#print('subsdatatsets', subDatasets)

ds = gdal.Warp('./output_lonlat.nc', src, **kwargs)
print("ds is ", ds)

del ds

Here is the command line operation that does work:

% gdalwarp -t_srs "+proj=longlat +datum=WGS84 +no_defs" -te -29 -2 -25 2 NETCDF:"./subset_20180729030000-JPL-L2P_GHRSST-SSTskin-MODIS_A-N-v02.0-fv01.0.nc":sea_surface_temperature -overwrite -of netCDF out_lonlat.nc
  • By the command line what you convert is literally this: NETCDF:"./subset_20180729030000-JPL-L2P_GHRSST-SSTskin-MODIS_A-N-v02.0-fv01.0.nc":sea_surface_temperature. The src that you use in the code is different. – user30184 Nov 13 '20 at 8:01
  • Yes the command line gdawarp only works on one netCDF layer (variable) at a time. A shortcoming in my mind. Cannot the python gdal.Warp operate on the entire netCDF file layer (variable) list all at once? Some examples I found seemed to indicate that. Perhaps I am mistaken, still new with this API. – Ed Armstrong Nov 13 '20 at 16:51
  • I do not believe that you can handle a bunch of subdatasets with gdalwarp. Gdal_translate does have a -sds option Copy all subdatasets of this file to individual output files. It should not be too difficult to do the same with Python: read the list of subdatasets and convert them one by one. – user30184 Nov 13 '20 at 17:51
  • You may be interested in: corteva.github.io/rioxarray/stable/examples/reproject.html – snowman2 Nov 15 '20 at 1:04
1

I found the solution:

  • I was using the arguments to outputBounds wrong. Its just an array and not a string.
  • As @user30184 suggested and I previously suspected gdal.Warp() only works on one netCDF layer (variable) at a time. To work on multiple variables you would need to read and execute on them one at a time.

Here is code that executes correctly:

        from osgeo import gdal, gdalconst
    
    
    kwargs = {'format': 'netCDF', 'copyMetadata': True, 
               'outputBounds': [-29, -2, -25, 2], 
              'dstSRS':'+proj=longlat +datum=WGS84 +no_defs',
             }
        
    # Source netCDF file
    src_filename = './subset_20180729030000-JPL-L2P_GHRSST-SSTskin-MODIS_A-N-v02.0-fv01.0.nc'
    
    # load the sea_surface_temperature variable
    nc_file = 'NETCDF:' + src_filename + ':sea_surface_temperature'
    src = gdal.Open(nc_file, gdalconst.GA_ReadOnly)
    
    # perform the warp and resample
    ds = gdal.Warp('./output_lonlat.nc', src, **kwargs)
    del ds
  • Another comment on gdal.Warp. Its does not completely copy the metadata from the source netCDF to the destination. Especially for the layer (variable) metadata, the often utilized scale_factor, add_offset, valid_min, valid_max are missing. These are critical to understand the layer and come up with the correct geophysical values. Annoying that they are not copied over. – Ed Armstrong Nov 15 '20 at 0:49
  • Another comment on gdal.Warp. Although it obviously understands CF coordinate systems and dimensions, the reference to CF as version 1.5 tells me the gdal maintainers have not looked at CF in a very long time. Current version is 1.8. Version 1.5 is from many years ago. – Ed Armstrong Nov 15 '20 at 0:51

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.