OGC WCS and TIME
Norman Barker
nbarker at ITTVIS.COM
Wed Dec 20 00:51:44 PST 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