OGC WCS and TIME

Norman Barker nbarker at ITTVIS.COM
Wed Dec 20 03:51:44 EST 2006


Hi Ara,

I have done this a lot :-), what type of datasets are you trying to use?
Python with GDAL, OGR makes this very easy.

For an example I have pasted a python script that creates U+V shape
indexes for MapServer WCS.  It prints the times out at the end which you
can copy in the mapfile as per the shape file. 

Let me know if you need any help.

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

	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:]

-----Original Message-----
From: UMN MapServer Users List [mailto:MAPSERVER-USERS at LISTS.UMN.EDU] On
Behalf Of Ara T Howard
Sent: 19 December 2006 23:30
To: MAPSERVER-USERS at LISTS.UMN.EDU
Subject: [UMN_MAPSERVER-USERS] OGC WCS and TIME

i'm trying to get time support running with wcs ala

   http://mapserver.gis.umn.edu/docs/howto/wms_time_support/

reading over

   http://mapserver.gis.umn.edu/docs/howto/wcs_server/#d45e385

and it seems like a little roll-your own shapefile munging code and and
appropriate metadata in the mapfile might do it.

has anyone does this already?  any examples out there?

regards.

-a
-- 
if you find yourself slandering anybody, first imagine that your mouth
is
filled with excrement.  it will break you of the habit quickly enough. -
the
dalai lama



More information about the mapserver-users mailing list