/[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 4044 - (hide annotations)
Tue Oct 30 03:35:17 2012 UTC (6 years, 10 months ago) by caltinay
File MIME type: text/x-python
File size: 34865 byte(s)
Some doco fixes for downunder.

1 caltinay 3946
2 jfenwick 3981 ##############################################################################
3 caltinay 3946 #
4     # Copyright (c) 2003-2012 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 caltinay 3946 #
7     # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10     #
11 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15 caltinay 3946
16 caltinay 4007 """Data readers/providers for inversions"""
17    
18 caltinay 3946 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
19 jfenwick 3981 http://www.uq.edu.au
20 caltinay 3946 Primary Business: Queensland, Australia"""
21     __license__="""Licensed under the Open Software License version 3.0
22     http://www.opensource.org/licenses/osl-3.0.php"""
23     __url__="https://launchpad.net/escript-finley"
24    
25 caltinay 3958 __all__ = ['DataSource','UBCDataSource','ERSDataSource','SyntheticDataSource','SmoothAnomaly']
26 caltinay 3947
27 caltinay 3946 import logging
28     import numpy as np
29 caltinay 4014 from esys.escript import ReducedFunction, Scalar
30 caltinay 4005 from esys.escript.linearPDEs import LinearSinglePDE
31 caltinay 4007 from esys.escript.util import *
32 caltinay 3946 import esys.escript.unitsSI as U
33 caltinay 4005 from esys.ripley import Brick, Rectangle, ripleycpp
34 caltinay 3971
35 caltinay 3947 try:
36     from scipy.io.netcdf import netcdf_file
37     __all__ += ['NetCDFDataSource']
38     except:
39     pass
40 caltinay 3946
41 caltinay 3957 def LatLonToUTM(lon, lat, wkt_string=None):
42     """
43     Converts one or more longitude,latitude pairs to the corresponding x,y
44     coordinates in the Universal Transverse Mercator projection.
45 caltinay 4007
46 caltinay 4020 :note: If the ``pyproj`` module is not installed a warning is printed and
47 caltinay 4007 the input values are scaled by a constant and returned.
48     :note: If `wkt_string` is not given or invalid or the ``gdal`` module is
49     not available to convert the string, then the input values are
50     assumed to be given in the Clarke 1866 projection.
51 caltinay 4009
52     :param lon: longitude value(s)
53     :type lon: `float`, `list`, `tuple`, or ``numpy.array``
54     :param lat: latitude value(s)
55     :type lat: `float`, `list`, `tuple`, or ``numpy.array``
56     :rtype: ``numpy.array``
57 caltinay 3957 """
58    
59     # not really optimal: if pyproj is not installed we return the input
60 caltinay 4005 # values scaled by a constant.
61 caltinay 3957 try:
62     import pyproj
63     except:
64     print("Warning, pyproj not available. Domain extents will be wrong")
65 caltinay 4009 return np.array(lon)*1000., np.array(lat)*1000.
66 caltinay 3957
67     # determine UTM zone from the input data
68 caltinay 3956 zone=int(np.median((np.floor((np.array(lon) + 180)/6) + 1) % 60))
69     try:
70     import osgeo.osr
71     srs = osgeo.osr.SpatialReference()
72     srs.ImportFromWkt(wkt_string)
73     p_src = pyproj.Proj(srs.ExportToProj4())
74     except:
75     p_src = pyproj.Proj('+proj=longlat +ellps=clrk66 +no_defs')
76 caltinay 3959 # we assume southern hemisphere here
77     p_dest = pyproj.Proj('+proj=utm +zone=%d +south +units=m'%zone)
78 caltinay 3956 x,y=pyproj.transform(p_src, p_dest, lon, lat)
79     return x,y
80 caltinay 3947
81 caltinay 3946 class DataSource(object):
82     """
83     A class that provides survey data for the inversion process.
84 caltinay 4005 This is an abstract base class that implements common functionality.
85     Methods to be overwritten by subclasses are marked as such.
86 caltinay 4044 This class assumes 2D data which is mapped to a slice of a 3D domain.
87     For other setups override the `_createDomain` method.
88 caltinay 3946 """
89    
90     def __init__(self):
91     """
92 caltinay 4005 Constructor. Sets some defaults and initializes logger.
93 caltinay 3946 """
94 caltinay 3958 self._constrainBottom=False
95 caltinay 3972 self._constrainSides=True
96 caltinay 3946 self._domain=None
97 gross 4039 self.__set_density_mask=None
98     self.__set_susceptibility_mask=None
99 caltinay 4005 self.setPadding()
100 caltinay 3946 self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
101    
102 caltinay 4005 def setPadding(self, pad_x=10, pad_y=10):
103 caltinay 3946 """
104 caltinay 4005 Sets the amount of padding around the dataset. If `pad_x`/`pad_y`
105     is >=1 the value is treated as number of elements to be added to the
106     domain (per side).
107     If ``0 < pad_x,pad_y < 1``, the padding amount is relative to the
108     dataset size. For example, calling ``setPadding(3, 0.1)`` to a data
109     source with size 10x20 will result in the padded data set size
110     16x24 (10+2*3, 20*(1+2*0.1))
111    
112     :param pad_x: Padding per side in x direction (default: 10 elements)
113     :type pad_x: ``int`` or ``float``
114     :param pad_y: Padding per side in y direction (default: 10 elements).
115     This value is only used for 3-dimensional datasets
116     :type pad_y: ``int`` or ``float``
117     """
118     self._pad_x=pad_x
119     self._pad_y=pad_y
120    
121     def setConstraints(self, bottom=False, sides=True):
122     """
123     If `bottom` is True, then the density mask will be set to 1 in the
124     padding area at the bottom of the domain. By default this area is
125     unconstrained. Similarly, if `sides` is True (default) then the
126     horizontal padding area is constrained, otherwise not.
127    
128     :param bottom: Whether to constrain the density at the bottom of the
129     domain
130     :type bottom: ``bool``
131     :param sides: Whether to constrain the density in the padding area
132     surrounding the data
133     :type sides: ``bool``
134     """
135     self._constrainBottom=bottom
136     self._constrainSides=sides
137    
138     def getDomain(self):
139     """
140     Returns a domain that spans the data area plus padding.
141     The domain is created the first time this method is called, subsequent
142     calls return the same domain so anything that affects the domain
143     (such as padding) needs to be set beforehand.
144    
145     :return: The escript domain for this data source.
146     :rtype: `esys.escript.Domain`
147     """
148     if self._domain is None:
149     self._domain=self._createDomain(self._pad_x, self._pad_y)
150     return self._domain
151    
152 gross 4032 def getSetDensityMask(self):
153 caltinay 4005 """
154     Returns the density mask data object, where mask has value 1 in the
155     padding area, 0 elsewhere.
156    
157     :return: The mask for the density.
158     :rtype: `esys.escript.Data`
159     """
160 gross 4039 return self.__set_density_mask
161 caltinay 4044
162 gross 4039 def setSetDensityMask(self, mask):
163 caltinay 4044 """
164     Sets the density mask to use.
165     """
166 gross 4039 self.__set_density_mask=mask
167 caltinay 4044
168 gross 4039 def getSetSusceptibilityMask(self):
169     """
170 caltinay 4044 Returns the susceptibility mask data object, where mask has value 1
171     in the padding area, 0 elsewhere.
172 caltinay 4005
173 gross 4039 :return: The mask for the susceptibility.
174     :rtype: `esys.escript.Data`
175     """
176     return self.__set_susceptibility_mask
177 caltinay 4044
178 gross 4039 def setSetSusceptibilityMask(self, mask):
179 caltinay 4044 """
180     Sets the susceptibility mask to use.
181     """
182 gross 4039 self.__set_susceptibility_mask=mask
183 caltinay 4044
184 caltinay 4005 def getGravityAndStdDev(self):
185     """
186     Returns the gravity anomaly and standard deviation data objects as a
187 caltinay 4044 tuple. This method must be implemented in subclasses that supply
188     gravity data.
189 caltinay 4005 """
190     raise NotImplementedError
191 caltinay 4044
192     def getMagneticFieldAndStdDev(self):
193 gross 4039 """
194 caltinay 4044 Returns the magnetic field and standard deviation data objects as a
195     tuple. This method must be implemented in subclasses that supply
196     magnetic data.
197 gross 4039 """
198     raise NotImplementedError
199 caltinay 4044
200 gross 4039 def getBackgroundMagneticField(self):
201     """
202 caltinay 4044 Returns the background magnetic field. This method must be
203     implemented in subclasses that supply magnetic data.
204 gross 4039 """
205     return NotImplementedError
206 caltinay 4005
207     def getDataExtents(self):
208     """
209     returns a tuple of tuples ``( (x0, y0), (nx, ny), (dx, dy) )``, where
210    
211     - ``x0``, ``y0`` = coordinates of data origin
212     - ``nx``, ``ny`` = number of data points in x and y
213     - ``dx``, ``dy`` = spacing of data points in x and y
214    
215     This method must be implemented in subclasses.
216     """
217     raise NotImplementedError
218    
219     def getVerticalExtents(self):
220     """
221     returns a tuple ``(z0, nz, dz)``, where
222    
223     - ``z0`` = minimum z coordinate (origin)
224 caltinay 4044 - ``nz`` = number of elements in z direction
225     - ``dz`` = spacing of elements (= cell size in z)
226 caltinay 4005
227     This method must be implemented in subclasses.
228     """
229     raise NotImplementedError
230    
231     def _addPadding(self, pad_x, pad_y, NE, l, origin):
232     """
233 caltinay 3946 Helper method that computes new number of elements, length and origin
234     after adding padding to the input values.
235 caltinay 4005
236     :param pad_x: Number of elements or fraction of padding in x direction
237     :type pad_x: ``int`` or ``float``
238     :param pad_y: Number of elements or fraction of padding in y direction
239     :type pad_y: ``int`` or ``float``
240     :param NE: Initial number of elements
241     :type NE: ``tuple`` or ``list``
242     :param l: Initial side lengths
243     :type l: ``tuple`` or ``list``
244     :param origin: Initial origin
245     :type origin: ``tuple`` or ``list``
246     :return: tuple with three elements ``(NE_padded, l_padded, origin_padded)``,
247     which are lists of the updated input parameters
248 caltinay 3946 """
249     DIM=len(NE)
250 caltinay 4005 frac=[0.]*(DIM-1)+[0]
251 caltinay 3946 # padding is applied to each side so multiply by 2 to get the total
252     # amount of padding per dimension
253 caltinay 4005 if pad_x>0 and pad_y<1:
254     frac[0]=2.*pad_x
255     elif pad_x>=1:
256     frac[0]=2.*pad_x/float(NE[0])
257     if DIM>2:
258     if pad_y>0 and pad_y<1:
259     frac[1]=2.*pad_y
260     elif pad_y>=1:
261     frac[1]=2.*pad_y/(float(NE[1]))
262 caltinay 3964
263 caltinay 3946 # calculate new number of elements
264 caltinay 4044 NE_new=[int(NE[i]*(1+frac[i])) for i in range(DIM)]
265     NEdiff=[NE_new[i]-NE[i] for i in range(DIM)]
266     spacing=[l[i]/NE[i] for i in range(DIM)]
267     l_new=[NE_new[i]*spacing[i] for i in range(DIM)]
268     origin_new=[origin[i]-NEdiff[i]/2.*spacing[i] for i in range(DIM)]
269 caltinay 3946 return NE_new, l_new, origin_new
270    
271 caltinay 3964 def _interpolateOnDomain(self, data):
272 caltinay 3946 """
273     Helper method that interpolates data arrays onto the domain.
274 caltinay 3965 Currently this works like a nearest neighbour mapping, i.e. values
275     are directly inserted into data objects at closest location.
276 caltinay 3946 """
277 caltinay 3964 dom=self.getDomain()
278     dim=dom.getDim()
279 caltinay 3965 # determine number of values required per element
280 caltinay 4014 DPP=Scalar(0., ReducedFunction(dom)).getNumberOfDataPoints()
281 caltinay 4044 for i in range(dim):
282 caltinay 3965 DPP=DPP/self._dom_NE[i]
283     DPP=int(DPP)
284    
285     # idx_mult.dot([x,y,z]) = flat index into data object
286     idx_mult=np.array([DPP]+self._dom_NE[:dim-1]).cumprod()
287    
288     # separate data arrays and coordinates
289     num_arrays=len(data[0])-dim
290     arrays=[]
291 caltinay 4044 for i in range(num_arrays):
292 caltinay 4014 d=Scalar(0., ReducedFunction(dom))
293 caltinay 3965 d.expand()
294     arrays.append(d)
295    
296     for entry in data:
297 caltinay 4044 index=[int((entry[i]-self._dom_origin[i])/self._spacing[i]) for i in range(dim)]
298 caltinay 3965 index=int(idx_mult.dot(index))
299 caltinay 4044 for i in range(num_arrays):
300     for p in range(DPP):
301 caltinay 3965 arrays[i].setValueOfDataPoint(index+p, entry[dim+i])
302    
303     return arrays
304    
305 caltinay 4007 def _createDomain(self, padding_x, padding_y):
306 caltinay 3965 """
307 caltinay 4005 Creates and returns an escript domain that spans the entire area of
308     available data plus a buffer zone. This method is called only once
309     the first time `getDomain()` is invoked and may be overwritten if
310     required.
311 caltinay 3946
312 caltinay 4005 :return: The escript domain for this data source.
313     :rtype: `esys.escript.Domain`
314 caltinay 3946 """
315 caltinay 3964 X0, NX, DX = self.getDataExtents()
316     z0, nz, dz = self.getVerticalExtents()
317 caltinay 3946
318 caltinay 3964 # number of elements (without padding)
319 caltinay 3965 NE = [NX[0], NX[1], nz]
320 caltinay 3946
321 caltinay 3964 # origin of domain (without padding)
322     origin = [X0[0], X0[1], z0]
323     origin = [np.round(oi) for oi in origin]
324 caltinay 3946
325 caltinay 3964 # cell size / point spacing
326     self._spacing = DX+[dz]
327 caltinay 3986 self._spacing = [float(np.round(si)) for si in self._spacing]
328 caltinay 3946
329 caltinay 3964 # length of domain (without padding)
330 caltinay 4044 l = [NE[i]*self._spacing[i] for i in range(len(NE))]
331 caltinay 3946
332 caltinay 3964 # now add padding to the values
333 caltinay 4007 NE_new, l_new, origin_new = self._addPadding(padding_x, padding_y, \
334 caltinay 3964 NE, l, origin)
335 caltinay 3946
336 caltinay 3964 # number of padding elements per side
337 caltinay 4044 NE_pad=[(NE_new[i]-NE[i])//2 for i in range(3)]
338 caltinay 3964
339 caltinay 3971 self._dom_NE_pad = NE_pad
340 caltinay 3964 self._dom_len = l_new
341     self._dom_NE = NE_new
342     self._dom_origin = origin_new
343 caltinay 4044 lo=[(origin_new[i], origin_new[i]+l_new[i]) for i in range(3)]
344     dom=Brick(*self._dom_NE, l0=lo[0], l1=lo[1], l2=lo[2])
345     # ripley may internally adjust NE and length, so recompute
346     self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(3)]
347     self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(3)]
348     x=dom.getX()-[self._dom_origin[i]+NE_pad[i]*self._spacing[i] for i in range(3)]
349     mask=wherePositive(dom.getX()[2])
350 caltinay 3946
351 caltinay 3964 # prepare density mask (=1 at padding area, 0 else)
352 caltinay 3972 if self._constrainSides:
353 caltinay 4044 for i in range(2):
354 caltinay 3972 mask=mask + whereNegative(x[i]) + \
355     wherePositive(x[i]-l_new[i]+2*NE_pad[i]*self._spacing[i])
356    
357 caltinay 3958 if self._constrainBottom:
358 caltinay 4009 mask = mask + whereNonPositive(x[2])
359 caltinay 4041 self.setSetDensityMask(wherePositive(mask))
360 caltinay 3946
361 caltinay 3964 self.logger.debug("Domain size: %d x %d x %d elements"%(self._dom_NE[0],self._dom_NE[1],self._dom_NE[2]))
362     self.logger.debug(" length: %g x %g x %g"%(self._dom_len[0],self._dom_len[1],self._dom_len[2]))
363 caltinay 3946 self.logger.debug(" origin: %g x %g x %g"%(origin_new[0],origin_new[1],origin_new[2]))
364    
365     return dom
366    
367 caltinay 3964
368     ##############################################################################
369     class UBCDataSource(DataSource):
370 caltinay 4044 def __init__(self, meshfile, gravfile, topofile=None):
371 caltinay 3964 super(UBCDataSource,self).__init__()
372     self.__meshfile=meshfile
373     self.__gravfile=gravfile
374     self.__topofile=topofile
375     self.__readMesh()
376    
377     def __readMesh(self):
378     meshdata=open(self.__meshfile).readlines()
379     numDataPoints=meshdata[0].split()
380     origin=meshdata[1].split()
381     self.__nPts=map(int, numDataPoints)
382     self.__origin=map(float, origin)
383     self.__delta=[float(X.split('*')[1]) for X in meshdata[2:]]
384     # vertical data is upside down
385     self.__origin[2]-=(self.__nPts[2]-1)*self.__delta[2]
386     self.logger.debug("Data Source: %s (mesh file: %s)"%(self.__gravfile, self.__meshfile))
387    
388     def getDataExtents(self):
389     """
390     returns ( (x0, y0), (nx, ny), (dx, dy) )
391     """
392     return (self.__origin[:2], self.__nPts[:2], self.__delta[:2])
393    
394     def getVerticalExtents(self):
395     """
396     returns (z0, nz, dz)
397     """
398     return (self.__origin[2], self.__nPts[2], self.__delta[2])
399    
400 gross 4032 #def getSetDensityMask(self):
401 caltinay 3964 # topodata=self.__readTopography()
402     # mask=self._interpolateOnDomain(topodata)
403     # mask=wherePositive(self.getDomain().getX()[2]-mask[0])
404     # return mask
405    
406     def getGravityAndStdDev(self):
407 caltinay 3966 gravlist=self.__readGravity() # x,y,z,g,s
408 caltinay 3964 g_and_sigma=self._interpolateOnDomain(gravlist)
409     return g_and_sigma[0]*[0,0,1], g_and_sigma[1]
410    
411     def __readTopography(self):
412     f=open(self.__topofile)
413 caltinay 3946 n=int(f.readline())
414     topodata=np.zeros((n,3))
415 caltinay 4044 for i in range(n):
416 caltinay 3946 x=f.readline().split()
417 caltinay 3964 x=map(float, x)
418 caltinay 3946 topodata[i]=x
419     f.close()
420     return topodata
421    
422 caltinay 3964 def __readGravity(self):
423     f=open(self.__gravfile)
424 caltinay 3946 n=int(f.readline())
425     gravdata=np.zeros((n,5))
426 caltinay 4044 for i in range(n):
427 caltinay 3946 x=f.readline().split()
428 caltinay 3964 x=map(float, x) # x, y, z, anomaly in mGal, stddev
429 caltinay 3946 # convert gravity anomaly units to m/s^2 and rescale error
430     x[3]*=-1e-5
431     x[4]*=1e-5
432     gravdata[i]=x
433     f.close()
434     return gravdata
435    
436     ##############################################################################
437     class NetCDFDataSource(DataSource):
438 caltinay 4044 """
439     Data Source for gridded netCDF data that use CF/COARDS conventions.
440     """
441 gross 4039 def __init__(self, gravfile=None, magfile=None, topofile=None, vertical_extents=(-40000,10000,25), alt_of_data=0.):
442 caltinay 3946 """
443 caltinay 3964 vertical_extents - (alt_min, alt_max, num_points)
444 caltinay 3946 alt_of_data - altitude of measurements
445     """
446     super(NetCDFDataSource,self).__init__()
447 caltinay 3964 self.__topofile=topofile
448     self.__gravfile=gravfile
449 gross 4039 self.__magfile=magfile
450 caltinay 3964 self.__determineExtents(vertical_extents)
451     self.__altOfData=alt_of_data
452 caltinay 3946
453 caltinay 3964 def __determineExtents(self, ve):
454     self.logger.debug("Data Source: %s"%self.__gravfile)
455     f=netcdf_file(self.__gravfile, 'r')
456 caltinay 3946 NX=0
457     for n in ['lon','longitude','x']:
458     if n in f.dimensions:
459     NX=f.dimensions[n]
460     break
461     if NX==0:
462     raise RuntimeError("Could not determine extents of data")
463     NY=0
464     for n in ['lat','latitude','y']:
465     if n in f.dimensions:
466     NY=f.dimensions[n]
467     break
468     if NY==0:
469     raise RuntimeError("Could not determine extents of data")
470    
471     # find longitude and latitude variables
472     lon_name=None
473     for n in ['lon','longitude']:
474 caltinay 3956 if n in f.variables:
475 caltinay 3946 lon_name=n
476 caltinay 3956 longitude=f.variables.pop(n)
477 caltinay 3946 break
478     if lon_name is None:
479     raise RuntimeError("Could not determine longitude variable")
480     lat_name=None
481     for n in ['lat','latitude']:
482 caltinay 3956 if n in f.variables:
483 caltinay 3946 lat_name=n
484 caltinay 3956 latitude=f.variables.pop(n)
485 caltinay 3946 break
486     if lat_name is None:
487     raise RuntimeError("Could not determine latitude variable")
488    
489     # try to figure out gravity variable name
490     grav_name=None
491 caltinay 3956 if len(f.variables)==1:
492     grav_name=f.variables.keys()[0]
493 caltinay 3946 else:
494 caltinay 3956 for n in f.variables.keys():
495     dims=f.variables[n].dimensions
496 caltinay 3946 if (lat_name in dims) and (lon_name in dims):
497     grav_name=n
498     break
499     if grav_name is None:
500     raise RuntimeError("Could not determine gravity variable")
501    
502 caltinay 4014 # try to determine value for unused data
503 caltinay 4016 if hasattr(f.variables[grav_name], 'missing_value'):
504 caltinay 4014 maskval = float(f.variables[grav_name].missing_value)
505 caltinay 4016 elif hasattr(f.variables[grav_name], '_FillValue'):
506     maskval = float(f.variables[grav_name]._FillValue)
507     else:
508     self.logger.debug("missing_value attribute not found, using default.")
509 caltinay 4014 maskval = 99999
510    
511 caltinay 3956 # see if there is a wkt string to convert coordinates
512     try:
513     wkt_string=f.variables[grav_name].esri_pe_string
514     except:
515     wkt_string=None
516    
517 caltinay 3958 # we don't trust actual_range & geospatial_lon_min/max since subset
518 caltinay 4014 # data does not seem to have these fields updated.
519     # Getting min/max from the arrays is obviously not very efficient but..
520 caltinay 3958 #lon_range=longitude.actual_range
521     #lat_range=latitude.actual_range
522     #lon_range=[f.geospatial_lon_min,f.geospatial_lon_max]
523     #lat_range=[f.geospatial_lat_min,f.geospatial_lat_max]
524     lon_range=longitude.data.min(),longitude.data.max()
525     lat_range=latitude.data.min(),latitude.data.max()
526     lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)
527     origin=[lon_range[0],lat_range[0],ve[0]]
528     lengths=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0],ve[1]-ve[0]]
529 caltinay 3956
530 caltinay 3946 f.close()
531 caltinay 3956
532 caltinay 3964 self.__nPts=[NX, NY, ve[2]]
533     self.__origin=origin
534 caltinay 4014 # we are rounding to avoid interpolation issues
535 caltinay 4044 self.__delta=[np.round(lengths[i]/self.__nPts[i]) for i in range(3)]
536 caltinay 3964 self.__wkt_string=wkt_string
537     self.__lon=lon_name
538     self.__lat=lat_name
539     self.__grv=grav_name
540 caltinay 4014 self.__maskval=maskval
541 caltinay 3946
542 caltinay 3964 def getDataExtents(self):
543     """
544     returns ( (x0, y0), (nx, ny), (dx, dy) )
545     """
546     return (self.__origin[:2], self.__nPts[:2], self.__delta[:2])
547 caltinay 3946
548 caltinay 3964 def getVerticalExtents(self):
549     """
550     returns (z0, nz, dz)
551     """
552     return (self.__origin[2], self.__nPts[2], self.__delta[2])
553 caltinay 3946
554 caltinay 3964 def getGravityAndStdDev(self):
555 caltinay 4014 nValues=self.__nPts[:2]+[1]
556     first=self._dom_NE_pad[:2]+[self._dom_NE_pad[2]+int((self.__altOfData-self.__origin[2])/self.__delta[2])]
557     g=ripleycpp._readNcGrid(self.__gravfile, self.__grv,
558     ReducedFunction(self.getDomain()),
559     first, nValues, (), self.__maskval)
560     sigma=whereNonZero(g-self.__maskval)
561     g=g*1e-6
562     sigma=sigma*2e-6
563     return g*[0,0,1], sigma
564 caltinay 3946
565     def _readTopography(self):
566 caltinay 3964 f=netcdf_file(self.__topofile, 'r')
567 caltinay 3946 lon=None
568     for n in ['lon','longitude']:
569     if n in f.variables:
570     lon=f.variables[n][:]
571     break
572     if lon is None:
573     raise RuntimeError("Could not determine longitude variable")
574     lat=None
575     for n in ['lat','latitude']:
576     if n in f.variables:
577     lat=f.variables[n][:]
578     break
579     if lat is None:
580     raise RuntimeError("Could not determine latitude variable")
581     alt=None
582     for n in ['altitude','alt']:
583     if n in f.variables:
584     alt=f.variables[n][:]
585     break
586     if alt is None:
587     raise RuntimeError("Could not determine altitude variable")
588    
589     topodata=np.column_stack((lon,lat,alt))
590     f.close()
591     return topodata
592    
593 caltinay 3958 ##############################################################################
594     class ERSDataSource(DataSource):
595     """
596     Data Source for ER Mapper raster data.
597     Note that this class only accepts a very specific type of ER Mapper data
598     input and will raise an exception if other data is found.
599     """
600 caltinay 4009 def __init__(self, headerfile, datafile=None, vertical_extents=(-40000,10000,25), alt_of_data=0.):
601 caltinay 3958 """
602     headerfile - usually ends in .ers
603     datafile - usually has the same name as the headerfile without '.ers'
604     """
605     super(ERSDataSource,self).__init__()
606 caltinay 3964 self.__headerfile=headerfile
607 caltinay 3958 if datafile is None:
608 caltinay 3964 self.__datafile=headerfile[:-4]
609 caltinay 3958 else:
610 caltinay 3964 self.__datafile=datafile
611     self.__readHeader(vertical_extents)
612     self.__altOfData=alt_of_data
613 caltinay 3958
614 caltinay 3964 def __readHeader(self, ve):
615     self.logger.debug("Data Source: %s (header: %s)"%(self.__datafile, self.__headerfile))
616     metadata=open(self.__headerfile, 'r').readlines()
617 caltinay 3958 # parse metadata
618     start=-1
619 caltinay 4044 for i in range(len(metadata)):
620 caltinay 3958 if metadata[i].strip() == 'DatasetHeader Begin':
621     start=i+1
622     if start==-1:
623     raise RuntimeError('Invalid ERS file (DatasetHeader not found)')
624    
625     md_dict={}
626     section=[]
627 caltinay 4044 for i in range(start, len(metadata)):
628 caltinay 3958 line=metadata[i]
629     if line[-6:].strip() == 'Begin':
630     section.append(line[:-6].strip())
631     elif line[-4:].strip() == 'End':
632     if len(section)>0:
633     section.pop()
634     else:
635     vals=line.split('=')
636     if len(vals)==2:
637     key = vals[0].strip()
638     value = vals[1].strip()
639     fullkey='.'.join(section+[key])
640     md_dict[fullkey]=value
641    
642     try:
643     if md_dict['RasterInfo.CellType'] != 'IEEE4ByteReal':
644     raise RuntimeError('Unsupported data type '+md_dict['RasterInfo.CellType'])
645     except KeyError:
646     self.logger.warn("Cell type not specified. Assuming IEEE4ByteReal.")
647    
648     try:
649     NX = int(md_dict['RasterInfo.NrOfCellsPerLine'])
650     NY = int(md_dict['RasterInfo.NrOfLines'])
651     except:
652     raise RuntimeError("Could not determine extents of data")
653    
654     try:
655     maskval = float(md_dict['RasterInfo.NullCellValue'])
656     except:
657 caltinay 3980 maskval = 99999
658 caltinay 3958
659     try:
660     spacingX = float(md_dict['RasterInfo.CellInfo.Xdimension'])
661     spacingY = float(md_dict['RasterInfo.CellInfo.Ydimension'])
662     except:
663     raise RuntimeError("Could not determine cell dimensions")
664    
665     try:
666 caltinay 3959 if md_dict['CoordinateSpace.CoordinateType']=='EN':
667     originX = float(md_dict['RasterInfo.RegistrationCoord.Eastings'])
668     originY = float(md_dict['RasterInfo.RegistrationCoord.Northings'])
669     elif md_dict['CoordinateSpace.CoordinateType']=='LL':
670     originX = float(md_dict['RasterInfo.RegistrationCoord.Longitude'])
671     originY = float(md_dict['RasterInfo.RegistrationCoord.Latitude'])
672     else:
673     raise RuntimeError("Unknown CoordinateType")
674 caltinay 3958 except:
675 caltinay 3959 self.logger.warn("Could not determine coordinate origin. Setting to (0.0, 0.0)")
676 caltinay 3958 originX,originY = 0.0, 0.0
677    
678 caltinay 3959 if 'GEODETIC' in md_dict['CoordinateSpace.Projection']:
679     # it appears we have lat/lon coordinates so need to convert
680     # origin and spacing. Try using gdal to get the wkt if available:
681     try:
682     from osgeo import gdal
683 caltinay 3964 ds=gdal.Open(self.__headerfile)
684 caltinay 3959 wkt=ds.GetProjection()
685     except:
686     wkt='GEOGCS["GEOCENTRIC DATUM of AUSTRALIA",DATUM["GDA94",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'
687 caltinay 3980 self.logger.warn('GDAL not available or file read error, assuming GDA94 data')
688 caltinay 3959 originX_UTM,originY_UTM=LatLonToUTM(originX, originY, wkt)
689     op1X,op1Y=LatLonToUTM(originX+spacingX, originY+spacingY, wkt)
690 caltinay 4009 # we are rounding to avoid interpolation issues
691 caltinay 3959 spacingX=np.round(op1X-originX_UTM)
692     spacingY=np.round(op1Y-originY_UTM)
693 caltinay 3964 originX=np.round(originX_UTM)
694     originY=np.round(originY_UTM)
695 caltinay 3959
696     # data sets have origin in top-left corner so y runs top-down
697 caltinay 3964 self.__dataorigin=[originX, originY]
698 caltinay 3958 originY-=(NY-1)*spacingY
699 caltinay 3964 self.__maskval=maskval
700 caltinay 4009 spacingZ=np.round(float(ve[1]-ve[0])/ve[2])
701 caltinay 3964 self.__delta = [spacingX, spacingY, spacingZ]
702     self.__nPts = [NX, NY, ve[2]]
703     self.__origin = [originX, originY, ve[0]]
704 caltinay 3958
705 caltinay 3964 def getDataExtents(self):
706     """
707     returns ( (x0, y0), (nx, ny), (dx, dy) )
708     """
709     return (self.__origin[:2], self.__nPts[:2], self.__delta[:2])
710 caltinay 3958
711 caltinay 3964 def getVerticalExtents(self):
712     """
713     returns (z0, nz, dz)
714     """
715     return (self.__origin[2], self.__nPts[2], self.__delta[2])
716    
717 caltinay 3958 def getGravityAndStdDev(self):
718 caltinay 3971 nValues=self.__nPts[:2]+[1]
719     first=self._dom_NE_pad[:2]+[self._dom_NE_pad[2]+int((self.__altOfData-self.__origin[2])/self.__delta[2])]
720 caltinay 3983 g=ripleycpp._readBinaryGrid(self.__datafile,
721     ReducedFunction(self.getDomain()),
722 caltinay 3980 first, nValues, (), self.__maskval)
723     sigma=whereNonZero(g-self.__maskval)
724 caltinay 3971 g=g*1e-6
725     sigma=sigma*2e-6
726     return g*[0,0,1], sigma
727 caltinay 3958
728    
729     ##############################################################################
730     class SourceFeature(object):
731     """
732     A feature adds a density distribution to (parts of) a domain of a synthetic
733     data source, for example a layer of a specific rock type or a simulated
734     ore body.
735     """
736     def getDensity(self):
737     """
738     Returns the density for the area covered by mask. It can be constant
739     or a data object with spatial dependency.
740     """
741     raise NotImplementedError
742     def getMask(self, x):
743     """
744     Returns the mask of the area of interest for this feature. That is,
745     mask is non-zero where the density returned by getDensity() should be
746     applied, zero elsewhere.
747     """
748     raise NotImplementedError
749    
750     class SmoothAnomaly(SourceFeature):
751 gross 4039 def __init__(self, lx, ly, lz, x, y, depth, rho_inner=None, rho_outer=None, k_inner=None, k_outer=None):
752 caltinay 3946 self.x=x
753     self.y=y
754     self.lx=lx
755     self.ly=ly
756     self.lz=lz
757     self.depth=depth
758     self.rho_inner=rho_inner
759     self.rho_outer=rho_outer
760 gross 4039 self.k_inner=k_inner
761     self.k_outer=k_outer
762 caltinay 3946 self.rho=None
763 gross 4039 self.k=None
764 caltinay 4044 self.mask=None
765 caltinay 4034
766 gross 4039 def getDensity(self,x):
767     if self.rho is None:
768 caltinay 4041 if self.rho_outer is None or self.rho_inner is None:
769     self.rho=0
770     else:
771 caltinay 4044 DIM=x.getDomain().getDim()
772 caltinay 4041 alpha=-log(abs(self.rho_outer/self.rho_inner))*4
773     rho=exp(-alpha*((x[0]-self.x)/self.lx)**2)
774     rho=rho*exp(-alpha*((x[DIM-1]-(sup(x[DIM-1])-self.depth))/self.lz)**2)
775     self.rho=maximum(abs(self.rho_outer), abs(self.rho_inner*rho))
776     if self.rho_inner<0: self.rho=-self.rho
777 caltinay 4044
778 caltinay 3958 return self.rho
779 caltinay 4044
780 gross 4039 def getSusceptibility(self,x):
781     if self.k is None:
782 caltinay 4041 if self.k_outer is None or self.k_inner is None:
783     self.k=0
784     else:
785 caltinay 4044 DIM=x.getDomain().getDim()
786 gross 4039 alpha=-log(abs(self.k_outer/self.k_inner))*4
787     k=exp(-alpha*((x[0]-self.x)/self.lx)**2)
788     k=k*exp(-alpha*((x[DIM-1]-(sup(x[DIM-1])-self.depth))/self.lz)**2)
789     self.k=maximum(abs(self.k_outer), abs(self.k_inner*k))
790     if self.k_inner<0: self.k=-self.k
791 caltinay 4044
792 gross 4039 return self.k
793 caltinay 4044
794 caltinay 3946 def getMask(self, x):
795     DIM=x.getDomain().getDim()
796     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)) \
797 caltinay 4034 *whereNonNegative(x[0]-(self.x-self.lx/2)) * whereNonPositive(x[0]-(self.x+self.lx/2))
798 caltinay 3946 if DIM>2:
799     m*=whereNonNegative(x[1]-(self.y-self.ly/2)) * whereNonPositive(x[1]-(self.y+self.ly/2))
800 caltinay 4044 self.mask = m
801 caltinay 4034 return m
802 caltinay 3946
803 caltinay 3958 ##############################################################################
804 caltinay 3946 class SyntheticDataSource(DataSource):
805 gross 4039 def __init__(self, DIM, NE, l, h, features, latitude=-32.):
806 caltinay 3946 super(SyntheticDataSource,self).__init__()
807     self._features = features
808     self.DIM=DIM
809     self.NE=NE
810     self.l=l
811     self.h=h
812 gross 4039 self.latitude=latitude
813 caltinay 3946
814 caltinay 4007 def _createDomain(self, padding_x, padding_y):
815 caltinay 3946 NE_H=self.NE
816     NE_L=int((self.l/self.h)*NE_H+0.5)
817     l=[self.l]*(self.DIM-1)+[self.h]
818     NE=[NE_L]*(self.DIM-1)+[NE_H]
819     origin=[0.]*self.DIM
820 caltinay 4007 NE_new, l_new, origin_new = self._addPadding(padding_x, padding_y, \
821 caltinay 3946 NE, l, origin)
822    
823     self.NE=NE_new
824     self.l=l_new[0]
825     self.h=l_new[self.DIM-1]
826    
827 caltinay 3960 self.logger.debug("Data Source: synthetic with %d features"%len(self._features))
828 caltinay 3946 if self.DIM==2:
829     dom = Rectangle(n0=NE_new[0], n1=NE_new[1], l0=l_new[0], l1=l_new[1])
830     self._x = dom.getX() + origin_new
831     self.logger.debug("Domain size: %d x %d elements"%(NE_new[0], NE_new[1]))
832     self.logger.debug(" length: %g x %g"%(l_new[0],l_new[1]))
833     self.logger.debug(" origin: %g x %g"%(origin_new[0],origin_new[1]))
834 caltinay 4034 else:
835 caltinay 3946 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])
836     self._x = dom.getX() + origin_new
837     self.logger.debug("Domain size: %d x %d x %d elements"%(self.NE[0],self.NE[1],self.NE[2]))
838     self.logger.debug(" length: %g x %g x %g"%(l_new[0],l_new[1],l_new[2]))
839     self.logger.debug(" origin: %g x %g x %g"%(origin_new[0],origin_new[1],origin_new[2]))
840    
841     dz=l_new[self.DIM-1]/NE_new[self.DIM-1]
842     self._g_mask=wherePositive(dom.getX()[0]-origin_new[0]) \
843     * whereNegative(dom.getX()[0]-(l_new[0]-origin_new[0])) \
844     * whereNonNegative(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+origin_new[self.DIM-1])) \
845     * whereNonPositive(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+(origin_new[self.DIM-1]+dz)))
846 caltinay 4044
847 gross 4039 self._B_mask=self._g_mask
848 caltinay 4044
849 gross 4039 mask=whereNegative(self._x[self.DIM-1]) + \
850 caltinay 3946 wherePositive(self._x[self.DIM-1]-l[self.DIM-1])
851 caltinay 4044 for i in range(self.DIM-1):
852 gross 4039 mask+= whereNegative(self._x[i]) + wherePositive(self._x[i]-l[i])
853     self.setSetDensityMask(wherePositive(mask))
854     self.setSetSusceptibilityMask(wherePositive(mask))
855 caltinay 3946
856     rho_ref=0.
857 gross 4039 k_ref=0
858 caltinay 3946 for f in self._features:
859     m=f.getMask(self._x)
860 gross 4039 rho_ref = rho_ref * (1-m) + f.getDensity(self._x) * m
861     k_ref = k_ref * (1-m) + f.getSusceptibility(self._x) * m
862 caltinay 3946 self._rho=rho_ref
863 gross 4039 self._k=k_ref
864 caltinay 3946
865     return dom
866    
867    
868     def getReferenceDensity(self):
869     return self._rho
870 gross 4039 def getReferenceSusceptibility(self):
871 caltinay 4044 return self._k
872 caltinay 3946
873     def getGravityAndStdDev(self):
874     pde=LinearSinglePDE(self.getDomain())
875 gross 4039 G=U.Gravitational_Constant
876 caltinay 3946 m_psi_ref=0.
877 caltinay 4044 for i in range(self.DIM):
878 caltinay 3946 m_psi_ref=m_psi_ref + whereZero(self._x[i]-inf(self._x[i])) \
879     + whereZero(self._x[i]-sup(self._x[i]))
880    
881     pde.setValue(A=kronecker(self.getDomain()), Y=-4*np.pi*G*self._rho, q=m_psi_ref)
882     pde.setSymmetryOn()
883     psi_ref=pde.getSolution()
884     del pde
885     g=-grad(psi_ref)
886     sigma=self._g_mask
887     return g,sigma
888 caltinay 4044
889     def getMagneticFieldAndStdDev(self):
890 gross 4039 pde=LinearSinglePDE(self.getDomain())
891     B_b=self.getBackgroundMagneticField()
892     DIM=self.getDomain().getDim()
893     m_psi_ref=0.
894 caltinay 4044 for i in range(self.DIM):
895 gross 4039 m_psi_ref=m_psi_ref + whereZero(self._x[i]-inf(self._x[i])) \
896 caltinay 4044 + whereZero(self._x[i]-sup(self._x[i]))
897 gross 4039 pde.setValue(A=kronecker(self.getDomain()), X=(1+self._k)*B_b, q=m_psi_ref)
898     pde.setSymmetryOn()
899     psi_ref=pde.getSolution()
900     del pde
901     B= (1+self._k) * B_b -grad(psi_ref)
902     sigma=self._B_mask
903     return B,sigma
904 caltinay 4044
905     def getBackgroundMagneticField(self):
906     theta = (90-self.latitude)/180.*np.pi
907     B_0=U.Mu_0 * U.Magnetic_Dipole_Moment_Earth / (4 * np.pi * U.R_Earth**3)
908     B_theta= B_0 * sin(theta)
909     B_r= 2 * B_0 * cos(theta)
910     DIM=self.getDomain().getDim()
911     if DIM<3:
912     return np.array([0., -B_r])
913     else:
914     return np.array([-B_theta, 0., -B_r])
915    

  ViewVC Help
Powered by ViewVC 1.1.26