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

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

fixed bounding box aggregation code

  • 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.minX=envelope.lowerCorner.CONTENT.split()[0]
51        self.minY=envelope.lowerCorner.CONTENT.split()[1]
52        self.maxX=envelope.upperCorner.CONTENT.split()[0]
53        self.maxY=envelope.upperCorner.CONTENT.split()[1]
54        self.t1= strToDate(envelope.beginPosition.CONTENT)
55        self.t2= strToDate(envelope.endPosition.CONTENT)
56 
57    def _compareLowerCorners(self,lowerCorner):
58        minX,minY=lowerCorner.CONTENT.split()[0],lowerCorner.CONTENT.split()[1]
59        if float(minX) < float(self.minX):
60            self.minX=minX
61        if float(minY) < float(self.minY):
62            self.minY=minY
63           
64    def _compareUpperCorners(self,upperCorner):
65        maxX,maxY=upperCorner.CONTENT.split()[0],upperCorner.CONTENT.split()[1]
66        if float(maxX) > float(self.maxX):
67            self.maxX=maxX
68        if float(maxY) > float(self.maxY):
69            self.maxY=maxY
70           
71    def _compareLowerTimes(self,timepos):
72        t=strToDate(timepos)
73        if t<self.t1:
74            self.t1=t
75       
76    def _compareUpperTimes(self,timepos):
77        t=strToDate(timepos)
78        if t>self.t2:
79            self.t2=t
80   
81    def compareEnvelope(self,envtocheck):
82        '''compares new envelope, and if necessary changes the aggregated envelope.'''
83        #if envtocheck.srsName!=self.envelope.srsName:
84            #print 'Skipping envelope with srs %s Can currently only perform aggregation within the same spatio-temporal reference system.'%envtocheck.srsName           
85        #else:
86        self._compareLowerCorners(envtocheck.lowerCorner)
87        self._compareUpperCorners(envtocheck.upperCorner)
88        self._compareLowerTimes(envtocheck.beginPosition.CONTENT)
89        self._compareUpperTimes(envtocheck.endPosition.CONTENT)
90   
91    def getAggregatedEnvelope(self):
92        newenvelope=csml.parser.EnvelopeWithTimePeriod()
93        newenvelope.beginPosition=csml.parser.csString(str(self.t1)) #convert from datetime types
94        newenvelope.endPosition=csml.parser.csString(str(self.t2)) 
95        newenvelope.lowerCorner=csml.parser.csString(str(self.minX +' '+ self.minY))
96        newenvelope.upperCorner=csml.parser.csString(str(self.maxX +' '+ self.maxY))
97        newenvelope.srsName='WGS84'
98        return newenvelope
99
100
101
102def getSeqRule(nDims):
103        #returns a sequenceRule
104        #*****************very simplified version*****************
105        #This really needs more work to see if the +s and -s are correct. and the x/y/z s in right order
106        #nDims = spatial dimensions + time dimension
107       
108        #TODO CSML2 - rewrite this
109        if nDims == 2:
110                seq = "Linear"
111        elif nDims ==3:
112                seq = "+x+y"
113        elif nDims ==4:
114                seq = "+x+y+z"
115        else:
116                seq="not sure"
117        return seq
118       
119       
120def getMappingRule(nDims):
121        #TODO CSML2 - rewrite this
122        if nDims == 2:
123                mr = "+ gridI + series"
124        elif nDims ==3:
125                mr = "+ gridI + gridJ +series"
126        elif nDims ==4:
127                mr = "+ gridI + gridJ +gridK +series"
128        else:
129                mr="not sure"
130        return mr
131       
132       
133def _getEnvelopeFromLimits(limits, axislabels):
134    env = csml.parser.EnvelopeWithTimePeriod()
135   
136    #simplistic reduction to remove time from CRS.
137    if limits['crs'][-1:] =='t':
138        env.srsName=limits['crs'][:-1]
139    else:
140        env.srsName=limits['crs']
141       
142    if env.srsName == 'ndg:crs:lonlat':
143        limits['crs'] ='WGS84'     
144        env.srsName = 'WGS84'     
145   
146    lowCorner= ''
147    upCorner=''
148    for axis in axislabels.split():
149        if axis != 'time':
150            try:
151                lowCorner = lowCorner + str(limits[axis][0]) + ' '
152                upCorner = upCorner + str(limits[axis][1]) + ' '
153            except:
154                pass
155        else:
156            env.beginPosition=csml.parser.csString(limits[axis][0])
157            env.endPosition=csml.parser.csString(limits[axis][1])
158    env.lowerCorner  = csml.parser.csString(lowCorner)
159    env.upperCorner  = csml.parser.csString(upCorner)
160    return env
161
162def addBoundingBoxes(ds):
163    '''Adds envelope with timeperiod to each feature, currently only works on GridSeriesFeatures but extendable.'''
164    timename='time'
165    cat=csml.csmllibs.csmlcrs.CRSCatalogue()
166    xlinkresolver=csml.csmllibs.csmlxlink.XlinkResolver(ds)
167    for feature in ds.featureCollection.featureMembers:
168        limits={}
169        if type(feature) is csml.parser.GridSeriesFeature:   
170            crsName=feature.value.gridSeriesDomain.srsName
171            limits['crs']=crsName
172            crs=cat.getCRS(feature.value.gridSeriesDomain.srsName)
173            axislabels=feature.value.gridSeriesDomain.axisLabels
174            for ord in feature.value.gridSeriesDomain.coordTransformTable.gridOrdinates:       
175                if ord.coordAxisLabel.CONTENT==timename:
176                    try:                       
177                        tmin= strToDate(ord.coordAxisValues.timePositionList.CONTENT.split()[0].split('T')[0])
178                        tmax=strToDate( ord.coordAxisValues.timePositionList.CONTENT.split()[len(ord.coordAxisValues.timePositionList.CONTENT.split())-1].split('T')[0])
179                    except AttributeError:
180                        sptlist=xlinkresolver.resolveXlink(ord.coordAxisValues.href[1:])                               
181                        try:
182                            tmin= strToDate(sptlist.timePositionList.CONTENT.split()[0].split('T')[0])
183                            tmax=strToDate(sptlist.timePositionList.CONTENT.split()[len(sptlist.timePositionList.CONTENT.split())-1].split('T')[0])
184                        except:
185                            tmin, tmax=0,0
186                    limits[timename]=[tmin,tmax]
187                else:
188                    if hasattr(ord.coordAxisValues, 'insertedExtract'):
189                        data=ord.coordAxisValues.insertedExtract.getData()[0].tolist()
190                    else:
191                        #TODO inline content
192                        #print ord.coordAxisValues.href
193                        fileextract= xlinkresolver.resolveXlink(ord.coordAxisValues.href[1:])                               
194                        data=fileextract.getData()[0].tolist()
195                    data.sort()
196                    #extend the bounding box to take account of cell bounds
197                    minval,maxval=_getMaxMinBounds(ord.coordAxisLabel.CONTENT,data)
198                    limits[ord.coordAxisLabel.CONTENT]=[minval,maxval]
199        feature.boundedBy=_getEnvelopeFromLimits(limits, axislabels)
200    return ds.featureCollection
201
202def _getMaxMinBounds(axisname,data):
203    if len(data) > 1:                   
204        halfcell=abs(0.5*(abs(data[1])-abs(data[0])))
205    else:
206        halfcell=0
207    minval=data[0] - halfcell
208    maxval=data[len(data)-1] + halfcell                   
209    if axisname=='longitude':
210        #deal with the case where the grid point is 0, so the grid goes from -x to +x
211        #note this only works if the grid is 0-360
212        if minval < 0:
213            if maxval - minval ==360.0:
214                minval=0.0
215                maxval=360.0
216    if axisname=='latitude':
217        #deal with latitudes
218        if minval < -90:
219            minval=-90.0
220        if maxval > 90:
221            maxval=90.0
222    return minval,maxval
223
224def addEnvelope(fc):
225    #adds EnvelopeWithTimePeriod to feature collection, based on the spatial extent of the features within.   
226    aggregator =  None
227    for feature in fc.featureMembers:
228        env=feature.boundedBy
229        if aggregator is None:
230            aggregator=EnvelopeAggregator(env)
231        else:
232            aggregator.compareEnvelope(env)
233    fcEnvelope=aggregator.getAggregatedEnvelope()
234    fc.boundedBy=fcEnvelope
235    return fc
236       
237               
238       
239# a couple of string cleaning functions:
240#   don't think these are used now
241def cleanString1(messystring):
242        #removes outer brackets and changes commas to spaces
243        cleanstring = messystring[1:-1]
244        cleanstring = cleanstring.replace(',',' ')
245        #strip off first (time) dimension.
246        #note, this might need rethinking for other datasets
247        cleanstring = cleanstring[3:]
248        return cleanstring
249#       
250def cleanString(messystring):
251   # removes outer brackets and 's, but leaves commas.
252    cleanstring = messystring[1:-1]
253    cleanstring = cleanstring.replace("'",'')
254    #remove any carriage returns
255    cleanstring = cleanstring.replace("\n",'')
256    return cleanstring
257       
258def listify(item):
259    ''' 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'''
260    if type(item) is list:
261        return item
262    else:
263        return [item]       
264
265def stringify(inputlist):
266    '''stringify takes a list and returns a string containing the contents of the list separated by spaces'''
267    strreturn=''
268    for item in inputlist:
269        strreturn=strreturn + item +' '
270    return strreturn[:-1]
271
272def checkDirExists(dirpath):
273    if os.path.isdir(dirpath) != True:
274        try:
275            os.mkdir(dirpath)
276        except:
277            pass
278             
279
280
281def main():
282    #test to check randomness of getRandomID()
283    idlist=[]
284    for j in range (500):
285        for i in range(100):
286            newid= getRandomID()
287            if newid in idlist:
288                print 'ID matched!!!'               
289                sys.exit()
290            else:               
291                idlist.append(newid)
292        print i * j
293    print 'done'
294   
295if __name__=="__main__":
296    main()
Note: See TracBrowser for help on using the repository browser.