/[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 3947 - (hide annotations)
Wed Aug 22 23:19:10 2012 UTC (7 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 20039 byte(s)
Compiling and installing downunder module now. Adjusted import statements
accordingly. Added a gravity test run.

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

  ViewVC Help
Powered by ViewVC 1.1.26