/[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 3946 - (hide annotations)
Wed Aug 22 02:08:12 2012 UTC (7 years ago) by caltinay
File MIME type: text/x-python
File size: 19896 byte(s)
Initial import of inversion code from devteam area, rev. 51.
Nothing is scons built/installed/tested yet.

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

  ViewVC Help
Powered by ViewVC 1.1.26