Loading NetCDF4 (HDF5) files into PostGIS Raster

A few months ago, I took a day or so to install version 2.0 of PostGIS and play around with it. In particular, I wanted to investigate the possibilities of storing our vast holdings of raster data (climate model output, hydrologic model output, gridded observations, and downscaled gridded observations) in the database, in the hopes of providing:

• a standard API for data selection
• spatial indexing on our rasters
• better ability to mix raster types and vector types (by use of all of PostGIS's geographic operators)
• etc., etc.

I got some data loaded, played around with viewing it in QGIS and then walked away from it for a couple months, without taking any notes or codifiying anything that I did. Bad idea. So today, I had to spend 30 minutes or so remembering why raster2pgsql (a data loader program which comes with PostGIS) couldn’t load my netcdfs. Weird.

james@basalt /tmp $ raster2pgsql -s 4326 -b 1-100 -a -f rast HDF5:"gfdl_cm2_1_20c_B1_daily.nc":tasmax bcsd | psql raster_test
ERROR: Unable to read raster file: HDF5:gfdl_cm2_1_20c_B1_daily.nc:tasmax

That’s weird, because gdalinfo can read it just fine; and raster2pgsql is supposedly just using GDAL…

james@basalt /tmp $ gdalinfo HDF5:"gfdl_cm2_1_20c_B1_daily.nc"://tasmax | head -100
Driver: HDF5Image/HDF5 Dataset
Files: gfdl_cm2_1_20c_B1_daily.nc
Size is 413, 199
Coordinate System is `'
Metadata:
  contact1=Markus Schnorbus
  contact2=mschnorb@uvic.ca
  contact3=tel: (250) 853-3502
  Conventions=CF-1.5
  experiment_id=gfdl_cm2.1_B1_run1
  grid_resolution=Square grid, 1/16 degrees_north by 1/16 degrees_east
  history=created 2011-10-19 13:41:34
  input_data=CMIP3 - GFDL CM2.1, B1 emissions scenario, run1
  institution=Pacific Climate Impacts Consortium (PCIC), Victoria, BC, www.pacificclimate.org
  project_id=Hydrologic Modelling
  reference=Werner (2011); Schnorbus et al. (2011)
  source=Bias Correction Spatial Disaggregation (BCSD) statistical downscaling
  spatial_domain=Irregular-shaped domain constrained to BC and northern edge of bordering US states
  title=Gridded meteorological data for British Columbia (BC), downscaled from GCM output
  version=1.0
  version_comment=derived from fortran symap implementation (non-PCIC)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  199.0)
Upper Right (  413.0,    0.0)
Lower Right (  413.0,  199.0)
Center      (  206.5,   99.5)
Band 1 Block=32x15 Type=Int16, ColorInterp=Undefined
  Metadata:
    tasmax_cell_method=time: maximum (interval: 1 day)
    tasmax_coordinates=lon lat
    tasmax_DIMENSION_LIST=
    tasmax_grid_mapping=geographic_wgs84
    tasmax_long_name=daily maximum surface air temperature
    tasmax_missing_value=-9999 
    tasmax_scale_factor=0.0099999998 
    tasmax_standard_name=surface_temperature
    tasmax_units=degrees_C
Band 2 Block=32x15 Type=Int16, ColorInterp=Undefined
  Metadata:
    tasmax_cell_method=time: maximum (interval: 1 day)
    tasmax_coordinates=lon lat
    tasmax_DIMENSION_LIST=
    tasmax_grid_mapping=geographic_wgs84
    tasmax_long_name=daily maximum surface air temperature
    tasmax_missing_value=-9999 
    tasmax_scale_factor=0.0099999998 
    tasmax_standard_name=surface_temperature
    tasmax_units=degrees_C

But wait… I remember this. google raster2pgsql hdf5 http://postgis.refractions.net/news/20120622/

Right… look at bugfix #1761. A bugfix for netcdfs/hdfs that was fixed in the first minor release. Sweet.

Unmask dev-db/postgis-2.0.1 in gentoo emerge -uv postgis

james@basalt /tmp $ raster2pgsql -s 4326 -b 1-100 -a -f rast HDF5:"./gfdl_cm2_1_20c_B1_daily.nc"://tasmax bcsd | psql raster_test
Processing 1/1: HDF5:./gfdl_cm2_1_20c_B1_daily.nc://tasmax

W00t.