/[escript]/trunk/downunder/py_src/datasources.py
ViewVC logotype

Annotation of /trunk/downunder/py_src/datasources.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3956 - (hide annotations)
Wed Sep 5 04:57:43 2012 UTC (7 years ago) by caltinay
File MIME type: text/x-python
File size: 20875 byte(s)
Fixes to netcdf reader for inversions.

1 caltinay 3946
2     ########################################################
3     #
4     # Copyright (c) 2003-2012 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7     #
8     # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11     #
12     ########################################################
13    
14     __copyright__="""Copyright (c) 2003-2012 by University of Queensland
15     Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18     __license__="""Licensed under the Open Software License version 3.0
19     http://www.opensource.org/licenses/osl-3.0.php"""
20     __url__="https://launchpad.net/escript-finley"
21    
22 caltinay 3947 __all__ = ['DataSource','UBCDataSource','SyntheticDataSource','SmoothAnomaly']
23    
24 caltinay 3946 import logging
25     import numpy as np
26 caltinay 3956 import pyproj
27 caltinay 3946 from esys.escript import *
28     from esys.escript.linearPDEs import *
29     import esys.escript.unitsSI as U
30 caltinay 3947 try:
31     from scipy.io.netcdf import netcdf_file
32     __all__ += ['NetCDFDataSource']
33     except:
34     pass
35 caltinay 3946
36 caltinay 3956 def LatLonToUTM(lon, lat, wkt_string):
37     zone=int(np.median((np.floor((np.array(lon) + 180)/6) + 1) % 60))
38     try:
39     import osgeo.osr
40     srs = osgeo.osr.SpatialReference()
41     srs.ImportFromWkt(wkt_string)
42     p_src = pyproj.Proj(srs.ExportToProj4())
43     except:
44     p_src = pyproj.Proj('+proj=longlat +ellps=clrk66 +no_defs')
45     p_dest = pyproj.Proj(proj='utm', zone=zone) # ellps?
46     x,y=pyproj.transform(p_src, p_dest, lon, lat)
47     return x,y
48 caltinay 3947
49 caltinay 3946 class DataSource(object):
50     """
51     A class that provides survey data for the inversion process.
52     """
53     # this is currently specific to gravity inversion and should be generalised
54    
55     def __init__(self):
56     """
57     """
58     self._domain=None
59     self._pad_l=0.1
60     self._pad_h=0.1
61     self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
62    
63     def _addPadding(self, pad_l, pad_h, NE, l, origin):
64     """
65     Helper method that computes new number of elements, length and origin
66     after adding padding to the input values.
67     """
68     DIM=len(NE)
69     frac=[0.]*DIM
70     # padding is applied to each side so multiply by 2 to get the total
71     # amount of padding per dimension
72     if pad_l>0 and pad_l<1:
73     for i in xrange(DIM-1):
74     frac[i]=2*pad_l
75     elif pad_l>=1:
76     for i in xrange(DIM-1):
77     frac[i]=2*pad_l/float(NE[i])
78     if pad_h>0 and pad_h<1:
79     frac[DIM-1]=2*pad_h
80     elif pad_h>=1:
81     frac[DIM-1]=2*pad_h/(float(NE[DIM-1]))
82     # calculate new number of elements
83     NE_new=[int(NE[i]*(1+frac[i])) for i in xrange(DIM)]
84     NEdiff=[NE_new[i]-NE[i] for i in xrange(DIM)]
85     spacing=[l[i]/NE[i] for i in xrange(DIM)]
86     l_new=[NE_new[i]*spacing[i] for i in xrange(DIM)]
87     origin_new=[origin[i]-NEdiff[i]/2.*spacing[i] for i in xrange(DIM)]
88     return NE_new, l_new, origin_new
89    
90     def _interpolateOnDomain(self, data, shape, origin, spacing, length):
91     """
92     Helper method that interpolates data arrays onto the domain.
93     """
94     dim=len(shape)
95     arrays=np.zeros(((len(data[0])-dim),)+tuple(shape))
96     for entry in data:
97     index=()
98     for i in range(dim):
99     index=((entry[i]-origin[i])/spacing[i],)+index
100     for i in range(arrays.shape[0]):
101     arrays[i][index]=entry[dim+i]
102     dom=self.getDomain()
103     x=dom.getX()
104     delta=[length[i]/(shape[dim-i-1]-1) for i in xrange(dim)]
105     realorigin=[inf(x[i]) for i in xrange(dim)]
106     res=[]
107     for i in range(arrays.shape[0]):
108     res.append(interpolateTable(arrays[i], x[:dim], realorigin, delta, 1e9))
109     return res
110    
111     def setPadding(self, pad_l=0.1, pad_h=0.1):
112     """
113     Sets the amount of padding around the dataset. If pad_l/pad_h is >=1
114     they are treated as number of elements to be added to the domain.
115     If 0 < pad_l;pad_h < 1, the padding amount is relative.
116     """
117     self._pad_l=pad_l
118     self._pad_h=pad_h
119    
120     def getDomain(self):
121     """
122     Returns a domain that spans the data area plus padding.
123     """
124     if self._domain is None:
125     self._domain=self._createDomain(self._pad_l, self._pad_h)
126     return self._domain
127    
128     def getDensityMask(self):
129     """
130     Returns the density mask data object, where mask has value 1 on the
131     padding area, 0 elsewhere.
132     """
133     raise NotImplementedError
134    
135     def getGravityAndStdDev(self):
136     """
137     Returns the gravity anomaly and standard deviation data objects as a
138     tuple.
139     """
140     raise NotImplementedError
141    
142     def _createDomain(self, padding_l, padding_h):
143     """
144     creates and returns an escript domain that spans the entire area of
145     available data plus a buffer zone.
146     """
147     raise NotImplementedError
148    
149    
150     ##############################################################################
151     class UBCDataSource(DataSource):
152     def __init__(self, domainclass, meshfile, gravfile, topofile=None):
153     super(UBCDataSource,self).__init__()
154     self._meshfile=meshfile
155     self._gravfile=gravfile
156     self._topofile=topofile
157     self._domainclass=domainclass
158     self._readMesh()
159    
160     def getDensityMask(self):
161     #topodata=self._readTopography()
162     #shape=[self.NE[1]+1, self.NE[0]+1]
163     #mask=self._interpolateOnDomain(topodata, shape, self._origin, self._spacing, self._meshlen)
164     #mask=wherePositive(self.getDomain().getX()[2]-mask[0])
165     return self._mask
166    
167     def getGravityAndStdDev(self):
168     gravlist=self._readGravity() # x,y,z,g,s
169     shape=[self.NE[2]+1, self.NE[1]+1, self.NE[0]+1]
170     g_and_sigma=self._interpolateOnDomain(gravlist, shape, self._origin, self._spacing, self._meshlen)
171     return g_and_sigma[0]*[0,0,1], g_and_sigma[1]
172    
173     def _readMesh(self):
174     meshdata=open(self._meshfile).readlines()
175     numDataPoints=meshdata[0].split()
176     origin=meshdata[1].split()
177     self._numDataPoints=[int(x) for x in numDataPoints]
178     self._origin=[float(x) for x in origin]
179     self._spacing=[float(X.split('*')[1]) for X in meshdata[2:]]
180     # vertical data is upside down
181     self._origin[2]-=(self._numDataPoints[2]-1)*self._spacing[2]
182    
183     def _createDomain(self, padding_l, padding_h):
184     NE=[self._numDataPoints[i]-1 for i in xrange(3)]
185     l=[NE[i]*self._spacing[i] for i in xrange(3)]
186     NE_new, l_new, origin_new = self._addPadding(padding_l, padding_h, \
187     NE, l, self._origin)
188    
189     self._meshlen=l_new
190     self.NE=NE_new
191     self._origin=origin_new
192     lo=[(self._origin[i], self._origin[i]+l_new[i]) for i in xrange(3)]
193     NEdiff=[NE_new[i]-NE[i] for i in xrange(3)]
194     try:
195     dom=self._domainclass(*self.NE, l0=lo[0], l1=lo[1], l2=lo[2])
196     x=dom.getX()-[self._origin[i]+NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]
197     mask=wherePositive(dom.getX()[2])
198    
199     except TypeError:
200     dom=self._domainclass(*self.NE, l0=l_new[0], l1=l_new[1], l2=l_new[2])
201     x=dom.getX()-[NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]
202     mask=wherePositive(x[2]+self._origin[2])
203    
204     M=2 # do not constrain bottom
205     #M=3 # constrain bottom
206     for i in xrange(M):
207     mask=mask + whereNegative(x[i]) + \
208     wherePositive(x[i]-l_new[i]+NEdiff[i]*self._spacing[i])
209     self._mask=wherePositive(mask)
210    
211     self.logger.debug("Domain size: %d x %d x %d elements"%(self.NE[0],self.NE[1],self.NE[2]))
212     self.logger.debug(" length: %g x %g x %g"%(l_new[0],l_new[1],l_new[2]))
213     self.logger.debug(" origin: %g x %g x %g"%(origin_new[0],origin_new[1],origin_new[2]))
214    
215     return dom
216    
217     def _readTopography(self):
218     f=open(self._topofile)
219     n=int(f.readline())
220     topodata=np.zeros((n,3))
221     for i in xrange(n):
222     x=f.readline().split()
223     x[0]=float(x[0])
224     x[1]=float(x[1])
225     x[2]=float(x[2])
226     topodata[i]=x
227     f.close()
228     return topodata
229    
230     def _readGravity(self):
231     f=open(self._gravfile)
232     n=int(f.readline())
233     gravdata=np.zeros((n,5))
234     for i in xrange(n):
235     x=f.readline().split()
236     x[0]=float(x[0]) # x
237     x[1]=float(x[1]) # y
238     x[2]=float(x[2]) # z
239     x[3]=float(x[3]) # gravity anomaly in mGal
240     x[4]=float(x[4]) # stddev
241     # convert gravity anomaly units to m/s^2 and rescale error
242     x[3]*=-1e-5
243     x[4]*=1e-5
244     gravdata[i]=x
245     f.close()
246     return gravdata
247    
248     ##############################################################################
249     class NetCDFDataSource(DataSource):
250 caltinay 3956 def __init__(self, domainclass, gravfile, topofile=None, vertical_extents=(-40000,10000,26), alt_of_data=1.):
251 caltinay 3946 """
252     vertical_extents - (alt_min, alt_max, num_elements)
253     alt_of_data - altitude of measurements
254     """
255     super(NetCDFDataSource,self).__init__()
256     self._topofile=topofile
257     self._gravfile=gravfile
258     self._domainclass=domainclass
259     self._determineExtents(vertical_extents)
260     self._altOfData=alt_of_data
261    
262     def getDensityMask(self):
263     #topodata=self._readTopography()
264     #shape=self._numDataPoints[1::-1]
265     #mask=self._interpolateOnDomain(topodata, shape, self._origin, self._spacing, self._meshlen)
266     #mask=wherePositive(self.getDomain().getX()[2]-mask[0])
267     #rho=fill*(1.-mask) + RHO_AIR*mask
268     return self._mask
269    
270     def getGravityAndStdDev(self):
271     gravlist=self._readGravity() # x,y,z,g,s
272     shape=[self.NE[2]+1, self.NE[1]+1, self.NE[0]+1]
273     g_and_sigma=self._interpolateOnDomain(gravlist, shape, self._origin, self._spacing, self._meshlen)
274     return g_and_sigma[0]*[0,0,1], g_and_sigma[1]
275    
276     def _determineExtents(self, ve):
277     f=netcdf_file(self._gravfile, 'r')
278     NX=0
279     for n in ['lon','longitude','x']:
280     if n in f.dimensions:
281     NX=f.dimensions[n]
282     break
283     if NX==0:
284     raise RuntimeError("Could not determine extents of data")
285     NY=0
286     for n in ['lat','latitude','y']:
287     if n in f.dimensions:
288     NY=f.dimensions[n]
289     break
290     if NY==0:
291     raise RuntimeError("Could not determine extents of data")
292    
293     # find longitude and latitude variables
294     lon_name=None
295     for n in ['lon','longitude']:
296 caltinay 3956 if n in f.variables:
297 caltinay 3946 lon_name=n
298 caltinay 3956 longitude=f.variables.pop(n)
299 caltinay 3946 break
300     if lon_name is None:
301     raise RuntimeError("Could not determine longitude variable")
302     lat_name=None
303     for n in ['lat','latitude']:
304 caltinay 3956 if n in f.variables:
305 caltinay 3946 lat_name=n
306 caltinay 3956 latitude=f.variables.pop(n)
307 caltinay 3946 break
308     if lat_name is None:
309     raise RuntimeError("Could not determine latitude variable")
310    
311     # try to figure out gravity variable name
312     grav_name=None
313 caltinay 3956 if len(f.variables)==1:
314     grav_name=f.variables.keys()[0]
315 caltinay 3946 else:
316 caltinay 3956 for n in f.variables.keys():
317     dims=f.variables[n].dimensions
318 caltinay 3946 if (lat_name in dims) and (lon_name in dims):
319     grav_name=n
320     break
321     if grav_name is None:
322     raise RuntimeError("Could not determine gravity variable")
323    
324 caltinay 3956 # see if there is a wkt string to convert coordinates
325     try:
326     wkt_string=f.variables[grav_name].esri_pe_string
327     except:
328     wkt_string=None
329    
330 caltinay 3946 origin=[0.,0.,ve[0]]
331 caltinay 3956 lengths=[100000.,100000.,ve[1]-ve[0]]
332    
333 caltinay 3946 try:
334 caltinay 3956 lon_range=longitude.actual_range
335     lat_range=latitude.actual_range
336     lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)
337     origin[:2]=lon_range[0],lat_range[0]
338 caltinay 3946 lengths[:2]=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]
339     except:
340     try:
341 caltinay 3956 lon_range=[f.geospatial_lon_min,f.geospatial_lon_max]
342     lat_range=[f.geospatial_lat_min,f.geospatial_lat_max]
343     lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)
344     origin[:2]=lon_range[0],lat_range[0]
345     lengths[:2]=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]
346 caltinay 3946 except:
347     pass
348    
349     f.close()
350 caltinay 3956
351 caltinay 3946 self._numDataPoints=[NX, NY, ve[2]]
352     self._origin=origin
353 caltinay 3956 self._spacing=[np.round(lengths[i]/(self._numDataPoints[i]-1)) for i in xrange(3)]
354     self._meshlen=[self._numDataPoints[i]*self._spacing[i] for i in xrange(3)]
355     self._wkt_string=wkt_string
356 caltinay 3946 self._lon=lon_name
357     self._lat=lat_name
358     self._grv=grav_name
359    
360     def _createDomain(self, padding_l, padding_h):
361     NE=[self._numDataPoints[i]-1 for i in xrange(3)]
362     l=self._meshlen
363     NE_new, l_new, origin_new = self._addPadding(padding_l, padding_h, \
364     NE, l, self._origin)
365    
366     self._meshlen=l_new
367     self.NE=NE_new
368     self._origin=origin_new
369     lo=[(self._origin[i], self._origin[i]+l_new[i]) for i in xrange(3)]
370     NEdiff=[NE_new[i]-NE[i] for i in xrange(3)]
371     try:
372     dom=self._domainclass(*self.NE, l0=lo[0], l1=lo[1], l2=lo[2])
373     # ripley may adjust NE and length
374     self._meshlen=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in xrange(3)]
375 caltinay 3956 self.NE=[self._meshlen[i]/self._spacing[i] for i in xrange(3)]
376 caltinay 3946 x=dom.getX()-[self._origin[i]+NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]
377     mask=wherePositive(dom.getX()[2])
378    
379     except TypeError:
380     dom=self._domainclass(*self.NE, l0=l_new[0], l1=l_new[1], l2=l_new[2])
381     x=dom.getX()-[NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]
382     mask=wherePositive(x[2]+self._origin[2])
383    
384     M=2 # do not constrain bottom
385     #M=3 # constrain bottom
386     for i in xrange(M):
387     mask=mask + whereNegative(x[i]) + \
388     wherePositive(x[i]-l_new[i]+NEdiff[i]*self._spacing[i])
389     self._mask=wherePositive(mask)
390    
391     self.logger.debug("Domain size: %d x %d x %d elements"%(self.NE[0],self.NE[1],self.NE[2]))
392     self.logger.debug(" length: %s x %s x %s"%(self._meshlen[0],self._meshlen[1],self._meshlen[2]))
393     self.logger.debug(" origin: %s x %s x %s"%(self._origin[0],self._origin[1],self._origin[2]))
394    
395     return dom
396    
397     def _readTopography(self):
398     f=netcdf_file(self._topofile, 'r')
399     lon=None
400     for n in ['lon','longitude']:
401     if n in f.variables:
402     lon=f.variables[n][:]
403     break
404     if lon is None:
405     raise RuntimeError("Could not determine longitude variable")
406     lat=None
407     for n in ['lat','latitude']:
408     if n in f.variables:
409     lat=f.variables[n][:]
410     break
411     if lat is None:
412     raise RuntimeError("Could not determine latitude variable")
413     alt=None
414     for n in ['altitude','alt']:
415     if n in f.variables:
416     alt=f.variables[n][:]
417     break
418     if alt is None:
419     raise RuntimeError("Could not determine altitude variable")
420    
421     topodata=np.column_stack((lon,lat,alt))
422     f.close()
423     return topodata
424    
425     def _readGravity(self):
426     f=netcdf_file(self._gravfile, 'r')
427 caltinay 3956 lon=f.variables[self._lon][:]
428     lat=f.variables[self._lat][:]
429     lon,lat=np.meshgrid(lon,lat)
430     lon,lat=LatLonToUTM(lon, lat, self._wkt_string)
431 caltinay 3946 grav=f.variables[self._grv][:]
432     f.close()
433     lon=lon.flatten()
434     lat=lat.flatten()
435     grav=grav.flatten()
436     alt=self._altOfData*np.ones(grav.shape)
437 caltinay 3956 # error value is an assumption
438 caltinay 3946 try:
439     missing=grav.missing_value
440 caltinay 3956 err=np.where(grav==missing, 20.0, 0.0)
441 caltinay 3946 except:
442 caltinay 3956 err=20.0*np.ones(lon.shape)
443     # convert units
444 caltinay 3946 err=1e-6*err
445     grav=1e-6*grav
446     gravdata=np.column_stack((lon,lat,alt,grav,err))
447     return gravdata
448    
449     class SmoothAnomaly(object):
450     def __init__(self, lx, ly, lz, x, y, depth, rho_inner, rho_outer):
451     self.x=x
452     self.y=y
453     self.lx=lx
454     self.ly=ly
455     self.lz=lz
456     self.depth=depth
457     self.rho_inner=rho_inner
458     self.rho_outer=rho_outer
459     self.rho=None
460     def getMask(self, x):
461     DIM=x.getDomain().getDim()
462     m=whereNonNegative(x[DIM-1]-(sup(x[DIM-1])-self.depth-self.lz/2)) * whereNonPositive(x[DIM-1]-(sup(x[DIM-1])-self.depth+self.lz/2)) \
463     *whereNonNegative(x[0]-(self.x-self.lx/2)) * whereNonPositive(x[0]-(self.x+self.lx/2))
464     if DIM>2:
465     m*=whereNonNegative(x[1]-(self.y-self.ly/2)) * whereNonPositive(x[1]-(self.y+self.ly/2))
466     if self.rho is None:
467     alpha=-log(abs(self.rho_outer/self.rho_inner))*4
468     rho=exp(-alpha*((x[0]-self.x)/self.lx)**2)
469     rho=rho*exp(-alpha*((x[DIM-1]-(sup(x[DIM-1])-self.depth))/self.lz)**2)
470     self.rho=maximum(abs(self.rho_outer), abs(self.rho_inner*rho))
471     if self.rho_inner<0: self.rho=-self.rho
472     return m
473    
474     class SyntheticDataSource(DataSource):
475     def __init__(self, DIM, NE, l, h, features):
476     super(SyntheticDataSource,self).__init__()
477     self._features = features
478     self.DIM=DIM
479     self.NE=NE
480     self.l=l
481     self.h=h
482    
483     def _createDomain(self, padding_l, padding_h):
484     NE_H=self.NE
485     NE_L=int((self.l/self.h)*NE_H+0.5)
486     l=[self.l]*(self.DIM-1)+[self.h]
487     NE=[NE_L]*(self.DIM-1)+[NE_H]
488     origin=[0.]*self.DIM
489     NE_new, l_new, origin_new = self._addPadding(padding_l, padding_h, \
490     NE, l, origin)
491    
492     self.NE=NE_new
493     self.l=l_new[0]
494     self.h=l_new[self.DIM-1]
495    
496     if self.DIM==2:
497     from esys.finley import Rectangle
498     dom = Rectangle(n0=NE_new[0], n1=NE_new[1], l0=l_new[0], l1=l_new[1])
499     self._x = dom.getX() + origin_new
500     self.logger.debug("Domain size: %d x %d elements"%(NE_new[0], NE_new[1]))
501     self.logger.debug(" length: %g x %g"%(l_new[0],l_new[1]))
502     self.logger.debug(" origin: %g x %g"%(origin_new[0],origin_new[1]))
503     else:
504     from esys.finley import Brick
505     dom = Brick(n0=NE_new[0], n1=NE_new[1], n2=NE_new[2], l0=l_new[0], l1=l_new[1], l2=l_new[2])
506     self._x = dom.getX() + origin_new
507     self.logger.debug("Domain size: %d x %d x %d elements"%(self.NE[0],self.NE[1],self.NE[2]))
508     self.logger.debug(" length: %g x %g x %g"%(l_new[0],l_new[1],l_new[2]))
509     self.logger.debug(" origin: %g x %g x %g"%(origin_new[0],origin_new[1],origin_new[2]))
510    
511     dz=l_new[self.DIM-1]/NE_new[self.DIM-1]
512     self._g_mask=wherePositive(dom.getX()[0]-origin_new[0]) \
513     * whereNegative(dom.getX()[0]-(l_new[0]-origin_new[0])) \
514     * whereNonNegative(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+origin_new[self.DIM-1])) \
515     * whereNonPositive(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+(origin_new[self.DIM-1]+dz)))
516     self._mask=whereNegative(self._x[self.DIM-1]) + \
517     wherePositive(self._x[self.DIM-1]-l[self.DIM-1])
518     for i in xrange(self.DIM-1):
519     self._mask=self._mask + whereNegative(self._x[i]) + \
520     wherePositive(self._x[i]-l[i])
521     self._mask=wherePositive(self._mask)
522    
523     rho_ref=0.
524     for f in self._features:
525     m=f.getMask(self._x)
526     rho_ref = rho_ref * (1-m) + f.rho * m
527     self._rho=rho_ref
528    
529     return dom
530    
531     def getDensityMask(self):
532     return self._mask
533    
534     def getReferenceDensity(self):
535     return self._rho
536    
537     def getGravityAndStdDev(self):
538     pde=LinearSinglePDE(self.getDomain())
539     G=6.6742e-11*U.m**3/(U.kg*U.sec**2)
540     m_psi_ref=0.
541     for i in range(self.DIM):
542     m_psi_ref=m_psi_ref + whereZero(self._x[i]-inf(self._x[i])) \
543     + whereZero(self._x[i]-sup(self._x[i]))
544    
545     pde.setValue(A=kronecker(self.getDomain()), Y=-4*np.pi*G*self._rho, q=m_psi_ref)
546     pde.setSymmetryOn()
547     psi_ref=pde.getSolution()
548     del pde
549     g=-grad(psi_ref)
550     sigma=self._g_mask
551     return g,sigma
552    

  ViewVC Help
Powered by ViewVC 1.1.26