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

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

Bounding boxes now correctly determined by scanner

  • 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                        fileextract= xlinkresolver.resolveXlink(ord.coordAxisValues.href[1:])                               
194                        data=fileextract.getData()[0].tolist()
195                    data.sort()
196                    minval=data[0]
197                    maxval=data[len(data)-1]
198                    limits[ord.coordAxisLabel.CONTENT]=[minval,maxval]
199        feature.boundedBy=_getEnvelopeFromLimits(limits, axislabels)
200    return ds.featureCollection
201
202
203
204
205
206def addEnvelope(fc, ffmap):
207    #adds EnvelopeWithTimePeriod to feature collection, based on the spatial extent of the features within.
208   
209    #process:   
210    #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.
211   
212    #TODO, need to use the CRS catalogue to get these values
213    lonname='longitude'
214    latname='latitude'
215    timename='time'
216    srs='WGS84'
217    tnewmin=None
218    tnewmax=None
219       
220    aggregator=None
221    for feature in fc.featureMembers:
222        if type(feature) is csml.parser.GridSeriesFeature:
223            tmax=None
224            tmin=None
225            for ord in feature.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
226                if ord.coordAxisLabel.CONTENT==timename:
227                    try:
228                        tnewmin= strToDate(ord.coordAxisValues.timePositionList.CONTENT.split()[0].split('T')[0])
229                        tnewmax=strToDate( ord.coordAxisValues.timePositionList.CONTENT.split()[len(ord.coordAxisValues.timePositionList.CONTENT.split())-1].split('T')[0])
230                    except AttributeError:
231                        pass  #using xlink, so these values will be picked up elsewhere.
232                    if tmin ==None:
233                        tmin=tnewmin
234                    if tmax == None:
235                        tmax=tnewmax
236                    if tnewmin < tmin:
237                        tmin=tnewmin                   
238                    if tnewmax > tmax:
239                        tmax=tnewmax                     
240        if type(feature) is csml.parser.ProfileSeriesFeature:
241            tmax=None
242            tmin=None
243            for ord in feature.value.profileSeriesDomain.coordTransformTable.gridOrdinates:
244                if ord.coordAxisLabel.CONTENT==timename:
245                    tnewmin= strToDate(ord.coordAxisValues.timePositionList.CONTENT.split()[0].split('T')[0])
246                    tnewmax=strToDate( ord.coordAxisValues.timePositionList.CONTENT.split()[len(ord.coordAxisValues.timePositionList.CONTENT.split())-1].split('T')[0])
247                    if tmin ==None:
248                        tmin=tnewmin
249                    if tmax == None:
250                        tmax=tnewmax
251                    if tnewmin < tmin:
252                        tmin=tnewmin                   
253                    if tnewmax > tmax:
254                        tmax=tnewmax                     
255                       
256    for repfile in ffmap.getRepresentativeFiles():
257        minlon=None
258        maxlon=None
259        minlat=None
260        maxlat=None
261        f=repfile.getRepresentativeFileName()
262        DI=csml.csmllibs.csmldataiface.DataInterface()
263        DI=DI.getUnknownInterfaceType(f)
264        DI.openFile(f)
265        tmpDims=DI.getListOfAxes()
266        if lonname in tmpDims:
267            DI.setAxis(lonname)
268            vals=DI.getDataForAxis()
269            minlon=str(min(vals))
270            maxlon= str(max(vals))
271            if latname in tmpDims:
272                DI.setAxis(latname)
273                vals=DI.getDataForAxis()
274                minlat=str(min(vals))
275                maxlat= str(max(vals))                         
276        else:
277            tmpVars=DI.getListofVariables()
278            if lonname in tmpVars:
279                DI.setVariable(lonname)
280                vals=DI.getDataForVar()
281                minlon=str(min(vals))
282                maxlon= str(max(vals))
283                if latname in tmpVars:
284                    DI.setVariable(latname)
285                    vals=DI.getDataForVar()
286                    minlat=str(min(vals))
287                    maxlat= str(max(vals))
288           
289        env=csml.parser.EnvelopeWithTimePeriod()
290        env.lowerCorner=csml.parser.csString(minlon + ' ' + minlat)
291        env.upperCorner=csml.parser.csString(maxlon + ' ' + maxlat)
292        env.beginPosition=csml.parser.csString(tmin)
293        env.endPosition= csml.parser.csString(tmax)
294        env.srsName=srs
295        if aggregator is None:
296            aggregator=EnvelopeAggregator(env)
297        else:
298            aggregator.compareEnvelope(env)
299    fcEnvelope=aggregator.getAggregatedEnvelope()
300    fcEnvelope.srsName=env.srsName 
301    fc.boundedBy=fcEnvelope
302    fc.boundedBy=env
303    return fc
304       
305       
306       
307       
308# a couple of string cleaning functions:
309#   don't think these are used now
310def cleanString1(messystring):
311        #removes outer brackets and changes commas to spaces
312        cleanstring = messystring[1:-1]
313        cleanstring = cleanstring.replace(',',' ')
314        #strip off first (time) dimension.
315        #note, this might need rethinking for other datasets
316        cleanstring = cleanstring[3:]
317        return cleanstring
318#       
319def cleanString(messystring):
320   # removes outer brackets and 's, but leaves commas.
321    cleanstring = messystring[1:-1]
322    cleanstring = cleanstring.replace("'",'')
323    #remove any carriage returns
324    cleanstring = cleanstring.replace("\n",'')
325    return cleanstring
326       
327def listify(item):
328    ''' 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'''
329    if type(item) is list:
330        return item
331    else:
332        return [item]       
333
334def stringify(inputlist):
335    '''stringify takes a list and returns a string containing the contents of the list separated by spaces'''
336    strreturn=''
337    for item in inputlist:
338        strreturn=strreturn + item +' '
339    return strreturn[:-1]
340
341def checkDirExists(dirpath):
342    if os.path.isdir(dirpath) != True:
343        try:
344            os.mkdir(dirpath)
345        except:
346            pass
347             
348
349
350def main():
351    #test to check randomness of getRandomID()
352    idlist=[]
353    for j in range (500):
354        for i in range(100):
355            newid= getRandomID()
356            if newid in idlist:
357                print 'ID matched!!!'               
358                sys.exit()
359            else:               
360                idlist.append(newid)
361        print i * j
362    print 'done'
363   
364if __name__=="__main__":
365    main()
Note: See TracBrowser for help on using the repository browser.