Simple python wrapper for space/time/variable subsetting with GDAL (geotiff) and NCO (netcdf)

[mpaget]

Define common parameters. Pass to GDAL and NCO depending on file type.

What output? Same as input type probably (different formats can be done by the reformatter)

[mgnethery]

Assorted thoughts from me.

The issue and possible solutions My initial problem to solve was to get a time-series for a single (latlon) point. I have also limited myself to thinking about two basic raster formats only (at this point): GeoTiff and netCDF. Various testing shows:

• for a netCDF file, ncks (from NCO package) works well, and also works well via DAP (ie. can query a Thredds/???? openDAP service). GDAL's 'gdallocationinfo' also works well.
• for a Thredds service, the NCSS (netCDF subsetting) request works well
• for a GeoTiff file, there is no easy way (I could find) to do a time series request. It's possible you could create a GDAL .vrt file and query using something like 'gdallocationinfo'.
• for GeoTiffs presented via GeoServer, the best option is to create an ImageMosaic layer. Further, to be able to query time-series info, I needed to install the 'netcdf-out" plugin. The request then becomes a WCS request, specifying netcdf as the output format.

So, considering the above and wanting to minimise confusion in having too many ways to query various formats, I have now concentrated on getting time-series via two methods:

• NCSS query for Thredds/netCDF
• WCS query for Geoserver/GeoTiff (with netcdf-out plugin)

Python script Based on this, I have written a python script ( ts-request.py ). It requires a yaml file ( auscover-products.yaml ) as a store of the products we want to expose. I envisage that a future cataloguer/inventoriser/crawler will be the store of this information. For starters though, the "ts-request.py" script (run interactively) does:

• display a list of products/layers we want to expose
• enter dimensions to subset on (spatial, temporal) and variable selection
• enter a geo-location
• extract time-series to a netCDF file, then reformat to csv

WPS A further issue to solve was to be able to make this query available via as web processing service (so click on a map, it will pass the geo-location to this WPS process and return a (csv) time-series for that point). Still early days on this but have implemented WPS via a plugin for GeoServer (plus a python scripting plugin). Intention is this specific WPS process will call a python broker, which will then call the aforementioned ts-request script. A simple broker script is located ( here ). An example WPS URL that is currently working: http://qld.auscover.org.au/geoserver-dev/ows?service=WPS&version=1.0.0&request=Execute&Identifier=py:timeSeries&dataInputs=layer=5;lon_lat_position=[geojson-string]

Summary of GeoServer

• most recent build
• plugins:
  • GDAL
  • JAI
  • netcdf-out
  • wps
  • python

[mgnethery]

Update: After much testing, Geoserver WPS plugin was too restrictive in not being able to specify metadata and default datatypes/mimetypes. Have now implemented the ZOO-project WPS service. This allows passing ComplexData, meaning geojson strings (eg. point location) can be passed to the WPS broker ( script ). A new example: http://qld.auscover.org.au/wps?Service=WPS&Request=Execute&Version=1.0.0&Identifier=timeSeries&dataInputs=layer=total_cover;lon_lat_position={%22type%22:%22FeatureCollection%22,%22features%22:[{%22type%22:%22Feature%22,%22geometry%22:{%22type%22:%22Point%22,%22coordinates%22:[149.00,-34.00,-2573.70]}}]}

[mpaget]

When selecting a grid from THREDDS' netcdfsubsetservice, it may fail with "unknown DataType == long". The reason is that somewhere in the backend ncss is choosing to write a netCDF3 files, which doesn't support the data type. The work-around would be to force a netCDF4 file to be written but not sure where this is would happen. Could submit an issue to https://github.com/Unidata/thredds.

The offending code will be somewhere near: thredds/tds/src/main/java/thredds/server/ncss/format The decision to write netcdf3 happens at or before: thredds/cdm/src/main/java/ucar/nc2/NetcdfFileWriter.java

[mgnethery]

Came across a very annoying issue which took ages to track down. It was specific to use of the rasterstats module (called by 'raster-stats.py').

Summary was I need to edit the file: /usr/lib/python2.7/ctypes/util.py", line 246, in _findSoname_ldconfig as per http://lists.osgeo.org/pipermail/zoo-discuss/2015-November/001382.html.