source: TI02-CSML/trunk/csml/csmllibs/csmlextra.py @ 3131

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/TI02-CSML/trunk/csml/csmllibs/csmlextra.py@3131
Revision 3131, 14.9 KB checked in by domlowe, 12 years ago (diff)

individual bounding boxes extended to edge of grid cells - not done yet for dataset bounding box

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1#Stuff with nowhere better to go
2import random, datetime
3import string
4import csml.parser
5import os, sys
6
7def getRandomID():
8    #returns a random 8 character alphanumeric string that can be used as an internal identifier in a CSML document.
9    #the ID only needs to be unique within the document that it occurs, so this is sufficiently random
10    #to make a clash extremely unlikely. The string is alphanumeric but always begins with an upper case letter. 
11    #Note: tested with 50000 ids and all were different
12    randomID=random.choice(string.letters)
13    randomID=string.upper(randomID)
14    for i in range(7):
15        randomID=randomID+random.choice(string.letters+string.digits)
16    return randomID
17
18def _locateByGmlId(elem, idtofind):
19    it=elem.getiterator()
20    for i in it:
21        if '{http://www.opengis.net/gml}id' in i.attrib:
22            if i.attrib['{http://www.opengis.net/gml}id'] == idtofind:
23                return i
24
25
26def _getEmptyCSMLobj(tag):
27    #given an elementtree tag name, such as SpatialOrTemporalPositionList
28    #return a parser object such as csml.parser.SpatialOrTemporalPosition()
29    csmlname=tag.split('}')[1]
30    return csml.parser.__dict__[csmlname]()       
31
32
33def getObjbyGMLID(elem,idtofind):
34    #Given an element tree object and a gml:id return the fully parsed CSML object
35    #that matches the id
36    #useful for finding xlink matches
37    matchingelement= _locateByGmlId(elem, idtofind)
38    csmlobj=_getEmptyCSMLobj(matchingelement.tag)
39    csmlobj.fromXML(matchingelement)
40    return csmlobj
41
42def strToDate(ymd):
43    '''given a string in form '2006-01-01' returns a datetime object'''
44    return datetime.date(*map(int, ymd.split('-')))
45
46class EnvelopeAggregator(object):
47    def __init__(self,envelope):
48        '''start with aggregated envelope equal to the initialising envelope
49        envelope must be of type csml.parser.EnvelopeWithTimePeriod'''
50        self.envelope=envelope
51        self.minX=self.envelope.lowerCorner.CONTENT.split()[0]
52        self.minY=self.envelope.lowerCorner.CONTENT.split()[1]
53        self.maxX=self.envelope.upperCorner.CONTENT.split()[0]
54        self.maxY=self.envelope.upperCorner.CONTENT.split()[1]
55        self.t1= strToDate(self.envelope.beginPosition.CONTENT)
56        self.t2= strToDate(self.envelope.endPosition.CONTENT)
57 
58    def _compareLowerCorners(self,lowerCorner):
59        minX,minY=lowerCorner.CONTENT.split()[0],lowerCorner.CONTENT.split()[1]
60        if float(minX) < float(self.minX):
61            self.envelope.lowerCorner.CONTENT=str(minX +' '+ self.minY)
62            self.minX=minX
63        if float(minY) < float(self.minY):
64            self.envelope.lowerCorner.CONTENT=str(self.minX +' '+ minY)
65            self.minY=minY
66           
67    def _compareUpperCorners(self,upperCorner):
68        maxX,maxY=upperCorner.CONTENT.split()[0],upperCorner.CONTENT.split()[1]
69        if float(maxX) > float(self.maxX):
70            self.envelope.upperCorner.CONTENT=str(maxX +' '+ self.maxY)
71            self.maxX=maxX
72        if float(maxY) > float(self.maxY):
73            self.envelope.upperCorner.CONTENT=str(self.maxX +' '+ maxY)
74            self.maxY=maxY
75           
76    def _compareLowerTimes(self,timepos):
77        t=strToDate(timepos)
78        if t<self.t1:
79            self.t1=t
80       
81    def _compareUpperTimes(self,timepos):
82        t=strToDate(timepos)
83        if t>self.t2:
84            self.t2=t
85   
86    def compareEnvelope(self,envtocheck):
87        '''compares new envelope, and if necessary changes the aggregated envelope.'''
88        if envtocheck.srsName!=self.envelope.srsName:
89            print 'Can currently only perform aggregation within the same spatio-temporal reference system.'
90            sys.exit()
91        else: 
92            self._compareLowerCorners(envtocheck.lowerCorner)
93            self._compareUpperCorners(envtocheck.upperCorner)
94            self._compareLowerTimes(envtocheck.beginPosition.CONTENT)
95            self._compareUpperTimes(envtocheck.endPosition.CONTENT)
96   
97    def getAggregatedEnvelope(self):
98        self.envelope.beginPosition=csml.parser.csString(str(self.t1)) #convert from datetime types
99        self.envelope.endPosition=csml.parser.csString(str(self.t2)) 
100        return self.envelope
101
102
103
104def getSeqRule(nDims):
105        #returns a sequenceRule
106        #*****************very simplified version*****************
107        #This really needs more work to see if the +s and -s are correct. and the x/y/z s in right order
108        #nDims = spatial dimensions + time dimension
109       
110        #TODO CSML2 - rewrite this
111        if nDims == 2:
112                seq = "Linear"
113        elif nDims ==3:
114                seq = "+x+y"
115        elif nDims ==4:
116                seq = "+x+y+z"
117        else:
118                seq="not sure"
119        return seq
120       
121       
122def getMappingRule(nDims):
123        #TODO CSML2 - rewrite this
124        if nDims == 2:
125                mr = "+ gridI + series"
126        elif nDims ==3:
127                mr = "+ gridI + gridJ +series"
128        elif nDims ==4:
129                mr = "+ gridI + gridJ +gridK +series"
130        else:
131                mr="not sure"
132        return mr
133       
134       
135def _getEnvelopeFromLimits(limits, axislabels):
136    env = csml.parser.EnvelopeWithTimePeriod()
137   
138    #simplistic reduction to remove time from CRS.
139    if limits['crs'][-1:] =='t':
140        env.srsName=limits['crs'][:-1]
141    else:
142        env.srsName=limits['crs']
143       
144    if env.srsName == 'ndg:crs:lonlat':
145        limits['crs'] ='WGS84'     
146        env.srsName = 'WGS84'     
147   
148    lowCorner= ''
149    upCorner=''
150    for axis in axislabels.split():
151        if axis != 'time':
152            try:
153                lowCorner = lowCorner + str(limits[axis][0]) + ' '
154                upCorner = upCorner + str(limits[axis][1]) + ' '
155            except:
156                pass
157        else:
158            env.beginPosition=csml.parser.csString(limits[axis][0])
159            env.endPosition=csml.parser.csString(limits[axis][1])
160    env.lowerCorner  = csml.parser.csString(lowCorner)
161    env.upperCorner  = csml.parser.csString(upCorner)
162    return env
163
164def addBoundingBoxes(ds):
165    '''Adds envelop with timeperiod to each feature, currently only works on GridSeriesFeatures but extendable.'''
166    timename='time'
167    cat=csml.csmllibs.csmlcrs.CRSCatalogue()
168    xlinkresolver=csml.csmllibs.csmlxlink.XlinkResolver(ds)
169    for feature in ds.featureCollection.featureMembers:
170        limits={}
171        if type(feature) is csml.parser.GridSeriesFeature:   
172            crsName=feature.value.gridSeriesDomain.srsName
173            limits['crs']=crsName
174            crs=cat.getCRS(feature.value.gridSeriesDomain.srsName)
175            axislabels=feature.value.gridSeriesDomain.axisLabels
176            for ord in feature.value.gridSeriesDomain.coordTransformTable.gridOrdinates:       
177                if ord.coordAxisLabel.CONTENT==timename:
178                    try:                       
179                        tmin= strToDate(ord.coordAxisValues.timePositionList.CONTENT.split()[0].split('T')[0])
180                        tmax=strToDate( ord.coordAxisValues.timePositionList.CONTENT.split()[len(ord.coordAxisValues.timePositionList.CONTENT.split())-1].split('T')[0])
181                    except AttributeError:
182                        sptlist=xlinkresolver.resolveXlink(ord.coordAxisValues.href[1:])                               
183                        try:
184                            tmin= strToDate(sptlist.timePositionList.CONTENT.split()[0].split('T')[0])
185                            tmax=strToDate(sptlist.timePositionList.CONTENT.split()[len(sptlist.timePositionList.CONTENT.split())-1].split('T')[0])
186                        except:
187                            tmin, tmax=0,0
188                    limits[timename]=[tmin,tmax]
189                else:
190                    if hasattr(ord.coordAxisValues, 'insertedExtract'):
191                        data=ord.coordAxisValues.insertedExtract.getData()[0].tolist()
192                    else:
193                        #TODO inline content
194                        #print ord.coordAxisValues.href
195                        fileextract= xlinkresolver.resolveXlink(ord.coordAxisValues.href[1:])                               
196                        data=fileextract.getData()[0].tolist()
197                    data.sort()
198                    #extend the bounding box to take account of cell bounds
199                    minval,maxval=_getMaxMinBounds(ord.coordAxisLabel.CONTENT,data)
200                    limits[ord.coordAxisLabel.CONTENT]=[minval,maxval]
201        feature.boundedBy=_getEnvelopeFromLimits(limits, axislabels)
202    return ds.featureCollection
203
204def _getMaxMinBounds(axisname,data):
205    if len(data) > 1:                   
206        halfcell=abs(0.5*(abs(data[1])-abs(data[0])))
207    else:
208        halfcell=0
209    minval=data[0] - halfcell
210    maxval=data[len(data)-1] + halfcell                   
211    if axisname=='longitude':
212        #deal with the case where the grid point is 0, so the grid goes from -x to +x
213        #note this only works if the grid is 0-360
214        if minval < 0:
215            if maxval - minval ==360.0:
216                minval=0.0
217                maxval=360.0
218    if axisname=='latitude':
219        #deal with latitudes
220        if minval < -90:
221            minval=-90.0
222        if maxval > 90:
223            maxval=90.0
224    return minval,maxval
225
226def addEnvelope(fc, ffmap):
227    #adds EnvelopeWithTimePeriod to feature collection, based on the spatial extent of the features within.
228   
229    #process:   
230    #Extract the lat/lon and time extents from each representative file.    Create a temporary EnvelopeWithTimePeriod object, then use the EnvelopeAggregator class to aggregate into an envelope for the whole feature collection.
231   
232    #TODO, need to use the CRS catalogue to get these values
233    lonname='longitude'
234    latname='latitude'
235    timename='time'
236    srs='WGS84'
237    tnewmin=None
238    tnewmax=None
239       
240    aggregator=None
241    for feature in fc.featureMembers:
242        if type(feature) is csml.parser.GridSeriesFeature:
243            tmax=None
244            tmin=None
245            for ord in feature.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
246                if ord.coordAxisLabel.CONTENT==timename:
247                    try:
248                        tnewmin= strToDate(ord.coordAxisValues.timePositionList.CONTENT.split()[0].split('T')[0])
249                        tnewmax=strToDate( ord.coordAxisValues.timePositionList.CONTENT.split()[len(ord.coordAxisValues.timePositionList.CONTENT.split())-1].split('T')[0])
250                    except AttributeError:
251                        pass  #using xlink, so these values will be picked up elsewhere.
252                    if tmin ==None:
253                        tmin=tnewmin
254                    if tmax == None:
255                        tmax=tnewmax
256                    if tnewmin < tmin:
257                        tmin=tnewmin                   
258                    if tnewmax > tmax:
259                        tmax=tnewmax                     
260        if type(feature) is csml.parser.ProfileSeriesFeature:
261            tmax=None
262            tmin=None
263            for ord in feature.value.profileSeriesDomain.coordTransformTable.gridOrdinates:
264                if ord.coordAxisLabel.CONTENT==timename:
265                    tnewmin= strToDate(ord.coordAxisValues.timePositionList.CONTENT.split()[0].split('T')[0])
266                    tnewmax=strToDate( ord.coordAxisValues.timePositionList.CONTENT.split()[len(ord.coordAxisValues.timePositionList.CONTENT.split())-1].split('T')[0])
267                    if tmin ==None:
268                        tmin=tnewmin
269                    if tmax == None:
270                        tmax=tnewmax
271                    if tnewmin < tmin:
272                        tmin=tnewmin                   
273                    if tnewmax > tmax:
274                        tmax=tnewmax                     
275                       
276    for repfile in ffmap.getRepresentativeFiles():
277        minlon=None
278        maxlon=None
279        minlat=None
280        maxlat=None
281        f=repfile.getRepresentativeFileName()
282        DI=csml.csmllibs.csmldataiface.DataInterface()
283        DI=DI.getUnknownInterfaceType(f)
284        DI.openFile(f)
285        tmpDims=DI.getListOfAxes()
286        if lonname in tmpDims:
287            DI.setAxis(lonname)
288            vals=DI.getDataForAxis()
289            minlon=str(min(vals))
290            maxlon= str(max(vals))
291            if latname in tmpDims:
292                DI.setAxis(latname)
293                vals=DI.getDataForAxis()
294                minlat=str(min(vals))
295                maxlat= str(max(vals))                         
296        else:
297            tmpVars=DI.getListofVariables()
298            if lonname in tmpVars:
299                DI.setVariable(lonname)
300                vals=DI.getDataForVar()
301                minlon=str(min(vals))
302                maxlon= str(max(vals))
303                if latname in tmpVars:
304                    DI.setVariable(latname)
305                    vals=DI.getDataForVar()
306                    minlat=str(min(vals))
307                    maxlat= str(max(vals))
308           
309        env=csml.parser.EnvelopeWithTimePeriod()
310        env.lowerCorner=csml.parser.csString(minlon + ' ' + minlat)
311        env.upperCorner=csml.parser.csString(maxlon + ' ' + maxlat)
312        env.beginPosition=csml.parser.csString(tmin)
313        env.endPosition= csml.parser.csString(tmax)
314        env.srsName=srs
315        if aggregator is None:
316            aggregator=EnvelopeAggregator(env)
317        else:
318            aggregator.compareEnvelope(env)
319    fcEnvelope=aggregator.getAggregatedEnvelope()
320    fcEnvelope.srsName=env.srsName 
321    fc.boundedBy=fcEnvelope
322    fc.boundedBy=env
323    return fc
324       
325       
326       
327       
328# a couple of string cleaning functions:
329#   don't think these are used now
330def cleanString1(messystring):
331        #removes outer brackets and changes commas to spaces
332        cleanstring = messystring[1:-1]
333        cleanstring = cleanstring.replace(',',' ')
334        #strip off first (time) dimension.
335        #note, this might need rethinking for other datasets
336        cleanstring = cleanstring[3:]
337        return cleanstring
338#       
339def cleanString(messystring):
340   # removes outer brackets and 's, but leaves commas.
341    cleanstring = messystring[1:-1]
342    cleanstring = cleanstring.replace("'",'')
343    #remove any carriage returns
344    cleanstring = cleanstring.replace("\n",'')
345    return cleanstring
346       
347def listify(item):
348    ''' listify checks if an item is a list, if it isn't it puts it inside a list and returns it. Always returns a list object'''
349    if type(item) is list:
350        return item
351    else:
352        return [item]       
353
354def stringify(inputlist):
355    '''stringify takes a list and returns a string containing the contents of the list separated by spaces'''
356    strreturn=''
357    for item in inputlist:
358        strreturn=strreturn + item +' '
359    return strreturn[:-1]
360
361def checkDirExists(dirpath):
362    if os.path.isdir(dirpath) != True:
363        try:
364            os.mkdir(dirpath)
365        except:
366            pass
367             
368
369
370def main():
371    #test to check randomness of getRandomID()
372    idlist=[]
373    for j in range (500):
374        for i in range(100):
375            newid= getRandomID()
376            if newid in idlist:
377                print 'ID matched!!!'               
378                sys.exit()
379            else:               
380                idlist.append(newid)
381        print i * j
382    print 'done'
383   
384if __name__=="__main__":
385    main()
Note: See TracBrowser for help on using the repository browser.