/[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 4014 - (hide annotations)
Thu Oct 4 03:28:35 2012 UTC (7 years ago) by caltinay
File MIME type: text/x-python
File size: 32144 byte(s)
Changed data source to use ripley's netCDF reader,
added tests and changed ERMapper test to have different x and y sizes.

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

  ViewVC Help
Powered by ViewVC 1.1.26