NetCDF WCS python script

[Norman Barker]

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:]
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/mapserver-users/attachments/20060728/99e1a77c/attachment.htm>