<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2//EN">
<HTML>
<HEAD>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=iso-8859-1">
<META NAME="Generator" CONTENT="MS Exchange Server version 6.5.7638.1">
<TITLE>NetCDF WCS python script</TITLE>
</HEAD>
<BODY>
<!-- Converted from text/plain format -->
<P><FONT SIZE=2>Hi,<BR>
<BR>
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.<BR>
<BR>
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.<BR>
<BR>
Norman<BR>
<BR>
# Create index shapefile<BR>
# derived from example at <A HREF="http://zcologia.com/news/categorylist_html?cat_id=2">http://zcologia.com/news/categorylist_html?cat_id=2</A> - Sean Gillies<BR>
import ogr<BR>
import os<BR>
import glob<BR>
import gdal<BR>
from time import strftime, strptime<BR>
<BR>
for var in ['u','v']:<BR>
<BR>
output = var + 'index.shp'<BR>
<BR>
# Create the Shapefile<BR>
driver = ogr.GetDriverByName('ESRI Shapefile')<BR>
if os.path.exists(output):<BR>
driver.DeleteDataSource(output)<BR>
<BR>
index_shp = driver.CreateDataSource(output)<BR>
index = index_shp.CreateLayer('quikscat',<BR>
geom_type=ogr.wkbPolygon)<BR>
<BR>
fd = ogr.FieldDefn('LOCATION', ogr.OFTString)<BR>
fd.SetWidth(100)<BR>
index.CreateField(fd)<BR>
<BR>
fdt = ogr.FieldDefn('TIME', ogr.OFTString)<BR>
fdt.SetWidth(30)<BR>
index.CreateField(fdt)<BR>
<BR>
# Loop over a number of georeferenced images<BR>
# The data files are 4 deep<BR>
files = glob.glob('*/*/*/*.nc')<BR>
times = ""<BR>
<BR>
for file in files:<BR>
# Get georeferencing and size of imagery<BR>
# read geoinfo from var component, missing geo-information at top level<BR>
subfilename = "NETCDF:" + "\"" + os.path.abspath(file) + "\"" + ":" + var;<BR>
dataset = gdal.Open(subfilename)<BR>
g = dataset.GetGeoTransform()<BR>
pixels = dataset.RasterXSize<BR>
lines = dataset.RasterYSize<BR>
minx = g[0]<BR>
maxx = minx + pixels * g[1]<BR>
maxy = g[3]<BR>
miny = maxy + lines * g[5]<BR>
<BR>
# append to the 'index' layer<BR>
wkt = 'POLYGON ((%f %f, %f %f, %f %f, %f %f, %f %f))' \<BR>
% (minx, miny, minx, maxy, maxx, maxy, maxx, miny, minx, miny)<BR>
g = ogr.CreateGeometryFromWkt(wkt)<BR>
<BR>
# calculate the date string, parse the ordinal date first<BR>
length = len(file)<BR>
dateStr = file[length - 12: length - 3] <BR>
date = strptime(dateStr, "%Y%j%H")<BR>
frmtDate = strftime("%Y-%m-%dT%H", date)<BR>
times += "," + frmtDate<BR>
f = ogr.Feature(feature_def=index.GetLayerDefn())<BR>
f.SetField(0, subfilename)<BR>
f.SetField(1, frmtDate)<BR>
f.SetGeometryDirectly(g)<BR>
index.CreateFeature(f)<BR>
f.Destroy()<BR>
<BR>
# destroy index_shp to flush and close<BR>
index_shp.Destroy()<BR>
<BR>
# print out times as they are useful for mapfile<BR>
print "TIMES VAR: " + var<BR>
print times[1:]<BR>
</FONT>
</P>
</BODY>
</HTML>