[UMN_MAPSERVER-USERS] NetCDF WCS python script

[Norman Barker]

Title: NetCDF WCS python script

Hi,

this isn't really suitable to go on the WCS how to (which uses Steve Lime's perl script as an example), but I have pasted a below a python script that creates mapserver time index files var1index.shp, var2index.shp etc., for the netCDF driver in GDAL with subdatasets.

The script is crude (it was hacked quickly), hopefully someone will find it useful.  The met community (as a result of GALEON) seem keen to host netCDF WCSs.

Norman

# Create index shapefile
# derived from example at http://zcologia.com/news/categorylist_html?cat_id=2 - Sean Gillies
import ogr
import os
import glob
import gdal
from time import strftime, strptime

for var in ['u','v']:

        output = var + 'index.shp'

        # Create the Shapefile
        driver = ogr.GetDriverByName('ESRI Shapefile')
        if os.path.exists(output):
                driver.DeleteDataSource(output)

        index_shp  = driver.CreateDataSource(output)
        index = index_shp.CreateLayer('quikscat',
                 geom_type=ogr.wkbPolygon)

        fd = ogr.FieldDefn('LOCATION', ogr.OFTString)
        fd.SetWidth(100)
        index.CreateField(fd)

        fdt = ogr.FieldDefn('TIME', ogr.OFTString)
        fdt.SetWidth(30)
        index.CreateField(fdt)

        # Loop over a number of georeferenced images
        # The data files are 4 deep
        files = glob.glob('*/*/*/*.nc')
        times = ""

        for file in files:
                # Get georeferencing and size of imagery
                # read geoinfo from var component, missing  geo-information at top level
                subfilename = "NETCDF:" + "\"" + os.path.abspath(file) + "\"" + ":" + var;
                dataset = gdal.Open(subfilename)
                g = dataset.GetGeoTransform()
                pixels = dataset.RasterXSize
                lines = dataset.RasterYSize
                minx = g[0]
                maxx = minx + pixels * g[1]
                maxy = g[3]
                miny = maxy + lines * g[5]
   
                # append to the 'index' layer
                wkt = 'POLYGON ((%f %f, %f %f, %f %f, %f %f, %f %f))' \
                        % (minx, miny, minx, maxy, maxx, maxy, maxx, miny, minx, miny)
                g = ogr.CreateGeometryFromWkt(wkt)

                # calculate the date string, parse the ordinal date first
                length = len(file)
                dateStr = file[length - 12: length - 3]   
                date = strptime(dateStr, "%Y%j%H")
                frmtDate = strftime("%Y-%m-%dT%H", date)
                times += "," + frmtDate
                f = ogr.Feature(feature_def=index.GetLayerDefn())
                f.SetField(0, subfilename)
                f.SetField(1, frmtDate)
                f.SetGeometryDirectly(g)
                index.CreateFeature(f)
                f.Destroy()

        # destroy index_shp to flush and close
        index_shp.Destroy()

        # print out times as they are useful for mapfile
        print "TIMES VAR: " + var
        print times[1:]