/[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 4016 - (hide annotations)
Tue Oct 9 03:50:27 2012 UTC (6 years, 11 months ago) by caltinay
File MIME type: text/x-python
File size: 32396 byte(s)
Skip data source tests under MPI for now as there is no straightforward way
of comparing the data (saveVTK appears to be the best option).

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 caltinay 4016 if hasattr(f.variables[grav_name], 'missing_value'):
479 caltinay 4014 maskval = float(f.variables[grav_name].missing_value)
480 caltinay 4016 elif hasattr(f.variables[grav_name], '_FillValue'):
481     maskval = float(f.variables[grav_name]._FillValue)
482     else:
483     self.logger.debug("missing_value attribute not found, using default.")
484 caltinay 4014 maskval = 99999
485    
486 caltinay 3956 # see if there is a wkt string to convert coordinates
487     try:
488     wkt_string=f.variables[grav_name].esri_pe_string
489     except:
490     wkt_string=None
491    
492 caltinay 3958 # we don't trust actual_range & geospatial_lon_min/max since subset
493 caltinay 4014 # data does not seem to have these fields updated.
494     # Getting min/max from the arrays is obviously not very efficient but..
495 caltinay 3958 #lon_range=longitude.actual_range
496     #lat_range=latitude.actual_range
497     #lon_range=[f.geospatial_lon_min,f.geospatial_lon_max]
498     #lat_range=[f.geospatial_lat_min,f.geospatial_lat_max]
499     lon_range=longitude.data.min(),longitude.data.max()
500     lat_range=latitude.data.min(),latitude.data.max()
501     lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)
502     origin=[lon_range[0],lat_range[0],ve[0]]
503     lengths=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0],ve[1]-ve[0]]
504 caltinay 3956
505 caltinay 3946 f.close()
506 caltinay 3956
507 caltinay 3964 self.__nPts=[NX, NY, ve[2]]
508     self.__origin=origin
509 caltinay 4014 # we are rounding to avoid interpolation issues
510     self.__delta=[np.round(lengths[i]/self.__nPts[i]) for i in xrange(3)]
511 caltinay 3964 self.__wkt_string=wkt_string
512     self.__lon=lon_name
513     self.__lat=lat_name
514     self.__grv=grav_name
515 caltinay 4014 self.__maskval=maskval
516 caltinay 3946
517 caltinay 3964 def getDataExtents(self):
518     """
519     returns ( (x0, y0), (nx, ny), (dx, dy) )
520     """
521     return (self.__origin[:2], self.__nPts[:2], self.__delta[:2])
522 caltinay 3946
523 caltinay 3964 def getVerticalExtents(self):
524     """
525     returns (z0, nz, dz)
526     """
527     return (self.__origin[2], self.__nPts[2], self.__delta[2])
528 caltinay 3946
529 caltinay 3964 def getDomainClass(self):
530     """
531     returns the domain generator class (e.g. esys.ripley.Brick)
532     """
533 caltinay 4014 return Brick
534 caltinay 3946
535 caltinay 3964 def getGravityAndStdDev(self):
536 caltinay 4014 nValues=self.__nPts[:2]+[1]
537     first=self._dom_NE_pad[:2]+[self._dom_NE_pad[2]+int((self.__altOfData-self.__origin[2])/self.__delta[2])]
538     g=ripleycpp._readNcGrid(self.__gravfile, self.__grv,
539     ReducedFunction(self.getDomain()),
540     first, nValues, (), self.__maskval)
541     sigma=whereNonZero(g-self.__maskval)
542     g=g*1e-6
543     sigma=sigma*2e-6
544     return g*[0,0,1], sigma
545 caltinay 3946
546     def _readTopography(self):
547 caltinay 3964 f=netcdf_file(self.__topofile, 'r')
548 caltinay 3946 lon=None
549     for n in ['lon','longitude']:
550     if n in f.variables:
551     lon=f.variables[n][:]
552     break
553     if lon is None:
554     raise RuntimeError("Could not determine longitude variable")
555     lat=None
556     for n in ['lat','latitude']:
557     if n in f.variables:
558     lat=f.variables[n][:]
559     break
560     if lat is None:
561     raise RuntimeError("Could not determine latitude variable")
562     alt=None
563     for n in ['altitude','alt']:
564     if n in f.variables:
565     alt=f.variables[n][:]
566     break
567     if alt is None:
568     raise RuntimeError("Could not determine altitude variable")
569    
570     topodata=np.column_stack((lon,lat,alt))
571     f.close()
572     return topodata
573    
574 caltinay 3958 ##############################################################################
575     class ERSDataSource(DataSource):
576     """
577     Data Source for ER Mapper raster data.
578     Note that this class only accepts a very specific type of ER Mapper data
579     input and will raise an exception if other data is found.
580     """
581 caltinay 4009 def __init__(self, headerfile, datafile=None, vertical_extents=(-40000,10000,25), alt_of_data=0.):
582 caltinay 3958 """
583     headerfile - usually ends in .ers
584     datafile - usually has the same name as the headerfile without '.ers'
585     """
586     super(ERSDataSource,self).__init__()
587 caltinay 3964 self.__headerfile=headerfile
588 caltinay 3958 if datafile is None:
589 caltinay 3964 self.__datafile=headerfile[:-4]
590 caltinay 3958 else:
591 caltinay 3964 self.__datafile=datafile
592     self.__readHeader(vertical_extents)
593     self.__altOfData=alt_of_data
594 caltinay 3958
595 caltinay 3964 def __readHeader(self, ve):
596     self.logger.debug("Data Source: %s (header: %s)"%(self.__datafile, self.__headerfile))
597     metadata=open(self.__headerfile, 'r').readlines()
598 caltinay 3958 # parse metadata
599     start=-1
600     for i in xrange(len(metadata)):
601     if metadata[i].strip() == 'DatasetHeader Begin':
602     start=i+1
603     if start==-1:
604     raise RuntimeError('Invalid ERS file (DatasetHeader not found)')
605    
606     md_dict={}
607     section=[]
608     for i in xrange(start, len(metadata)):
609     line=metadata[i]
610     if line[-6:].strip() == 'Begin':
611     section.append(line[:-6].strip())
612     elif line[-4:].strip() == 'End':
613     if len(section)>0:
614     section.pop()
615     else:
616     vals=line.split('=')
617     if len(vals)==2:
618     key = vals[0].strip()
619     value = vals[1].strip()
620     fullkey='.'.join(section+[key])
621     md_dict[fullkey]=value
622    
623     try:
624     if md_dict['RasterInfo.CellType'] != 'IEEE4ByteReal':
625     raise RuntimeError('Unsupported data type '+md_dict['RasterInfo.CellType'])
626     except KeyError:
627     self.logger.warn("Cell type not specified. Assuming IEEE4ByteReal.")
628    
629     try:
630     NX = int(md_dict['RasterInfo.NrOfCellsPerLine'])
631     NY = int(md_dict['RasterInfo.NrOfLines'])
632     except:
633     raise RuntimeError("Could not determine extents of data")
634    
635     try:
636     maskval = float(md_dict['RasterInfo.NullCellValue'])
637     except:
638 caltinay 3980 maskval = 99999
639 caltinay 3958
640     try:
641     spacingX = float(md_dict['RasterInfo.CellInfo.Xdimension'])
642     spacingY = float(md_dict['RasterInfo.CellInfo.Ydimension'])
643     except:
644     raise RuntimeError("Could not determine cell dimensions")
645    
646     try:
647 caltinay 3959 if md_dict['CoordinateSpace.CoordinateType']=='EN':
648     originX = float(md_dict['RasterInfo.RegistrationCoord.Eastings'])
649     originY = float(md_dict['RasterInfo.RegistrationCoord.Northings'])
650     elif md_dict['CoordinateSpace.CoordinateType']=='LL':
651     originX = float(md_dict['RasterInfo.RegistrationCoord.Longitude'])
652     originY = float(md_dict['RasterInfo.RegistrationCoord.Latitude'])
653     else:
654     raise RuntimeError("Unknown CoordinateType")
655 caltinay 3958 except:
656 caltinay 3959 self.logger.warn("Could not determine coordinate origin. Setting to (0.0, 0.0)")
657 caltinay 3958 originX,originY = 0.0, 0.0
658    
659 caltinay 3959 if 'GEODETIC' in md_dict['CoordinateSpace.Projection']:
660     # it appears we have lat/lon coordinates so need to convert
661     # origin and spacing. Try using gdal to get the wkt if available:
662     try:
663     from osgeo import gdal
664 caltinay 3964 ds=gdal.Open(self.__headerfile)
665 caltinay 3959 wkt=ds.GetProjection()
666     except:
667     wkt='GEOGCS["GEOCENTRIC DATUM of AUSTRALIA",DATUM["GDA94",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'
668 caltinay 3980 self.logger.warn('GDAL not available or file read error, assuming GDA94 data')
669 caltinay 3959 originX_UTM,originY_UTM=LatLonToUTM(originX, originY, wkt)
670     op1X,op1Y=LatLonToUTM(originX+spacingX, originY+spacingY, wkt)
671 caltinay 4009 # we are rounding to avoid interpolation issues
672 caltinay 3959 spacingX=np.round(op1X-originX_UTM)
673     spacingY=np.round(op1Y-originY_UTM)
674 caltinay 3964 originX=np.round(originX_UTM)
675     originY=np.round(originY_UTM)
676 caltinay 3959
677     # data sets have origin in top-left corner so y runs top-down
678 caltinay 3964 self.__dataorigin=[originX, originY]
679 caltinay 3958 originY-=(NY-1)*spacingY
680 caltinay 3964 self.__maskval=maskval
681 caltinay 4009 spacingZ=np.round(float(ve[1]-ve[0])/ve[2])
682 caltinay 3964 self.__delta = [spacingX, spacingY, spacingZ]
683     self.__nPts = [NX, NY, ve[2]]
684     self.__origin = [originX, originY, ve[0]]
685 caltinay 3958
686 caltinay 3964 def getDataExtents(self):
687     """
688     returns ( (x0, y0), (nx, ny), (dx, dy) )
689     """
690     return (self.__origin[:2], self.__nPts[:2], self.__delta[:2])
691 caltinay 3958
692 caltinay 3964 def getVerticalExtents(self):
693     """
694     returns (z0, nz, dz)
695     """
696     return (self.__origin[2], self.__nPts[2], self.__delta[2])
697    
698     def getDomainClass(self):
699     """
700     returns the domain generator class (e.g. esys.ripley.Brick)
701     """
702 caltinay 3971 return Brick
703 caltinay 3964
704 caltinay 3958 def getGravityAndStdDev(self):
705 caltinay 3971 nValues=self.__nPts[:2]+[1]
706     first=self._dom_NE_pad[:2]+[self._dom_NE_pad[2]+int((self.__altOfData-self.__origin[2])/self.__delta[2])]
707 caltinay 3983 g=ripleycpp._readBinaryGrid(self.__datafile,
708     ReducedFunction(self.getDomain()),
709 caltinay 3980 first, nValues, (), self.__maskval)
710     sigma=whereNonZero(g-self.__maskval)
711 caltinay 3971 g=g*1e-6
712     sigma=sigma*2e-6
713     return g*[0,0,1], sigma
714 caltinay 3958
715    
716     ##############################################################################
717     class SourceFeature(object):
718     """
719     A feature adds a density distribution to (parts of) a domain of a synthetic
720     data source, for example a layer of a specific rock type or a simulated
721     ore body.
722     """
723     def getDensity(self):
724     """
725     Returns the density for the area covered by mask. It can be constant
726     or a data object with spatial dependency.
727     """
728     raise NotImplementedError
729     def getMask(self, x):
730     """
731     Returns the mask of the area of interest for this feature. That is,
732     mask is non-zero where the density returned by getDensity() should be
733     applied, zero elsewhere.
734     """
735     raise NotImplementedError
736    
737     class SmoothAnomaly(SourceFeature):
738 caltinay 3946 def __init__(self, lx, ly, lz, x, y, depth, rho_inner, rho_outer):
739     self.x=x
740     self.y=y
741     self.lx=lx
742     self.ly=ly
743     self.lz=lz
744     self.depth=depth
745     self.rho_inner=rho_inner
746     self.rho_outer=rho_outer
747     self.rho=None
748 caltinay 3958
749     def getDensity(self):
750     return self.rho
751    
752 caltinay 3946 def getMask(self, x):
753     DIM=x.getDomain().getDim()
754     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)) \
755     *whereNonNegative(x[0]-(self.x-self.lx/2)) * whereNonPositive(x[0]-(self.x+self.lx/2))
756     if DIM>2:
757     m*=whereNonNegative(x[1]-(self.y-self.ly/2)) * whereNonPositive(x[1]-(self.y+self.ly/2))
758     if self.rho is None:
759     alpha=-log(abs(self.rho_outer/self.rho_inner))*4
760     rho=exp(-alpha*((x[0]-self.x)/self.lx)**2)
761     rho=rho*exp(-alpha*((x[DIM-1]-(sup(x[DIM-1])-self.depth))/self.lz)**2)
762     self.rho=maximum(abs(self.rho_outer), abs(self.rho_inner*rho))
763     if self.rho_inner<0: self.rho=-self.rho
764     return m
765    
766 caltinay 3958 ##############################################################################
767 caltinay 3946 class SyntheticDataSource(DataSource):
768     def __init__(self, DIM, NE, l, h, features):
769     super(SyntheticDataSource,self).__init__()
770     self._features = features
771     self.DIM=DIM
772     self.NE=NE
773     self.l=l
774     self.h=h
775    
776 caltinay 4007 def _createDomain(self, padding_x, padding_y):
777 caltinay 3946 NE_H=self.NE
778     NE_L=int((self.l/self.h)*NE_H+0.5)
779     l=[self.l]*(self.DIM-1)+[self.h]
780     NE=[NE_L]*(self.DIM-1)+[NE_H]
781     origin=[0.]*self.DIM
782 caltinay 4007 NE_new, l_new, origin_new = self._addPadding(padding_x, padding_y, \
783 caltinay 3946 NE, l, origin)
784    
785     self.NE=NE_new
786     self.l=l_new[0]
787     self.h=l_new[self.DIM-1]
788    
789 caltinay 3960 self.logger.debug("Data Source: synthetic with %d features"%len(self._features))
790 caltinay 3946 if self.DIM==2:
791     from esys.finley import Rectangle
792     dom = Rectangle(n0=NE_new[0], n1=NE_new[1], l0=l_new[0], l1=l_new[1])
793     self._x = dom.getX() + origin_new
794     self.logger.debug("Domain size: %d x %d elements"%(NE_new[0], NE_new[1]))
795     self.logger.debug(" length: %g x %g"%(l_new[0],l_new[1]))
796     self.logger.debug(" origin: %g x %g"%(origin_new[0],origin_new[1]))
797     else:
798     from esys.finley import Brick
799     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])
800     self._x = dom.getX() + origin_new
801     self.logger.debug("Domain size: %d x %d x %d elements"%(self.NE[0],self.NE[1],self.NE[2]))
802     self.logger.debug(" length: %g x %g x %g"%(l_new[0],l_new[1],l_new[2]))
803     self.logger.debug(" origin: %g x %g x %g"%(origin_new[0],origin_new[1],origin_new[2]))
804    
805     dz=l_new[self.DIM-1]/NE_new[self.DIM-1]
806     self._g_mask=wherePositive(dom.getX()[0]-origin_new[0]) \
807     * whereNegative(dom.getX()[0]-(l_new[0]-origin_new[0])) \
808     * whereNonNegative(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+origin_new[self.DIM-1])) \
809     * whereNonPositive(dom.getX()[self.DIM-1]-(l_new[self.DIM-1]+(origin_new[self.DIM-1]+dz)))
810     self._mask=whereNegative(self._x[self.DIM-1]) + \
811     wherePositive(self._x[self.DIM-1]-l[self.DIM-1])
812     for i in xrange(self.DIM-1):
813     self._mask=self._mask + whereNegative(self._x[i]) + \
814     wherePositive(self._x[i]-l[i])
815     self._mask=wherePositive(self._mask)
816    
817     rho_ref=0.
818     for f in self._features:
819     m=f.getMask(self._x)
820 caltinay 3958 rho_ref = rho_ref * (1-m) + f.getDensity() * m
821 caltinay 3946 self._rho=rho_ref
822    
823     return dom
824    
825     def getDensityMask(self):
826     return self._mask
827    
828     def getReferenceDensity(self):
829     return self._rho
830    
831     def getGravityAndStdDev(self):
832     pde=LinearSinglePDE(self.getDomain())
833     G=6.6742e-11*U.m**3/(U.kg*U.sec**2)
834     m_psi_ref=0.
835 caltinay 3964 for i in xrange(self.DIM):
836 caltinay 3946 m_psi_ref=m_psi_ref + whereZero(self._x[i]-inf(self._x[i])) \
837     + whereZero(self._x[i]-sup(self._x[i]))
838    
839     pde.setValue(A=kronecker(self.getDomain()), Y=-4*np.pi*G*self._rho, q=m_psi_ref)
840     pde.setSymmetryOn()
841     psi_ref=pde.getSolution()
842     del pde
843     g=-grad(psi_ref)
844     sigma=self._g_mask
845     return g,sigma
846    

  ViewVC Help
Powered by ViewVC 1.1.26