/[escript]/trunk/downunder/py_src/domainbuilder.py
ViewVC logotype

Annotation of /trunk/downunder/py_src/domainbuilder.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4613 - (hide annotations)
Mon Jan 13 01:09:01 2014 UTC (5 years, 9 months ago) by caltinay
File MIME type: text/x-python
File size: 19822 byte(s)
Addresses #43 - altitude is now scaled according to reference system height
unit.
Also separated out UTM zone determination from the actual projection code
to silence bogus output.

1 caltinay 4060
2     ##############################################################################
3     #
4 caltinay 4145 # Copyright (c) 2003-2013 by University of Queensland
5 caltinay 4060 # http://www.uq.edu.au
6     #
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     # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15    
16     """Domain construction from survey data for inversions"""
17    
18 caltinay 4145 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
19 caltinay 4060 http://www.uq.edu.au
20     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     __all__ = ['DomainBuilder']
26    
27     import logging
28     import numpy as np
29     from esys.escript.util import *
30 caltinay 4097 from esys.escript import unitsSI as U
31 caltinay 4060 from esys.ripley import Brick, Rectangle
32     from .datasources import DataSource
33 gross 4373 from .coordinates import ReferenceSystem, CartesianReferenceSystem
34 caltinay 4060
35     class DomainBuilder(object):
36     """
37     This class is responsible for constructing an escript Domain object with
38     suitable extents and resolution for survey data (`DataSource` objects)
39 caltinay 4435 that are added to it.
40    
41     The domain covers a region above and below the Earth surface. The
42     East-West direction is used as the x- or longitudinal or x[0] direction,
43     the North-South direction is used as the y- or latitudinal or x[1]
44     direction, the vertical direction is denoted by z or radial or x[2]
45     direction. The corresponding terms are used synonymously.
46 caltinay 4060 """
47    
48 gross 4373 def __init__(self, dim=3, reference_system=None):
49 caltinay 4060 """
50     Constructor.
51    
52     :param dim: Dimensionality (2 or 3) of the target domain.
53     This has implications for the survey data than can be
54     added. By default a 3D domain is created.
55     :type dim: ``int``
56 gross 4373 :param reference_system: reference coordinate system. By default the
57     Cartesian coordinate system is used.
58     :type reference_system: `ReferenceSystem`
59 caltinay 4060 """
60     self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
61     if dim not in (2,3):
62     raise ValueError("Number of dimensions must be 2 or 3")
63 caltinay 4394 if not reference_system:
64     self.__reference_system=CartesianReferenceSystem()
65     else:
66     self.__reference_system=reference_system
67 caltinay 4435
68 caltinay 4436 if self.__reference_system.isCartesian():
69     self.__v_scale=1.
70     else:
71     self.__v_scale=1./self.getReferenceSystem().getHeightUnit()
72    
73 caltinay 4145 self.__domain=None
74     self.__dim=dim
75     self.__sources=[]
76     self.__background_magnetic_field=None
77 gross 4373 self.setElementPadding()
78 caltinay 4060 self.setVerticalExtents()
79 gross 4120 self.fixDensityBelow()
80     self.fixSusceptibilityBelow()
81 gross 4373
82     def getReferenceSystem(self):
83     """
84     returns the reference coordinate system
85    
86     :rtype: `ReferenceSystem`
87     """
88     return self.__reference_system
89    
90 caltinay 4060 def addSource(self, source):
91     """
92     Adds a survey data provider to the domain builder.
93 caltinay 4201 An exception is raised if the domain has already been built or if the
94     UTM zone of `source` does not match the UTM zone of sources already
95 caltinay 4362 added to the domain builder (see Inversion Cookbook for more
96     information). An exception is also raised if the dimensionality of the
97     data source is incompatible with this domain builder. That is, the
98     dimensionality of the data must be one less than the dimensionality
99     of the domain (specified in the constructor).
100 caltinay 4060
101 gross 4373 :param source: The data source to be added. Its reference system needs
102     to match the reference system of the DomainBuilder.
103 caltinay 4060 :type source: `DataSource`
104     """
105 caltinay 4150 if self.__domain is not None:
106     raise RuntimeError("Invalid call to addSource(). Domain is already built.")
107 caltinay 4060 if not isinstance(source, DataSource):
108     raise TypeError("source is not a DataSource")
109 caltinay 4394 if not source.getReferenceSystem() == self.getReferenceSystem():
110     raise ValueError("source reference system does not match.")
111 caltinay 4201
112 caltinay 4362 DATA_DIM = len(source.getDataExtents()[0])
113     if DATA_DIM != self.__dim-1:
114     raise ValueError("Data must be %d-dimensional."%(self.__dim-1))
115 caltinay 4201 if len(self.__sources)>0:
116     if self.__sources[0].getUtmZone() != source.getUtmZone():
117     raise ValueError("It is not possible to combine data sources located in different UTM zones at the moment.")
118    
119 caltinay 4145 self.__sources.append(source)
120 caltinay 4060
121 gross 4373 def setFractionalPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
122 caltinay 4060 """
123 caltinay 4108 Sets the amount of padding around the dataset as a fraction of the
124     dataset side lengths.
125    
126 caltinay 4150 For example, calling ``setFractionalPadding(0.2, 0.1)`` with a data
127     source of size 10x20 will result in the padded data set size
128 caltinay 4108 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
129 caltinay 4060
130 caltinay 4108 :param pad_x: Padding per side in x direction (default: no padding)
131     :type pad_x: ``float``
132     :param pad_y: Padding per side in y direction (default: no padding)
133     :type pad_y: ``float``
134 gross 4373 :param pad_lat: Padding per side in latitudinal direction (default: no padding)
135     :type pad_lat: ``float``
136     :param pad_lon: Padding per side in longitudinal direction (default: no padding)
137     :type pad_lon: ``float``
138     :note: `pad_y` is ignored for 2-dimensional domains.
139 caltinay 4060 """
140 gross 4373 if not pad_lat == None:
141 caltinay 4394 if not pad_x == None:
142     raise ValueError("Either pad_lat or pad_x can be set.")
143     else:
144     pad_x = pad_lat
145 gross 4373 if not pad_lon == None:
146 caltinay 4394 if not pad_y == None:
147     raise ValueError("Either pad_lon or pad_y can be set.")
148     else:
149     pad_y = pad_lan
150 caltinay 4150 if self.__domain is not None:
151     raise RuntimeError("Invalid call to setFractionalPadding(). Domain is already built.")
152 caltinay 4108 if pad_x is not None:
153     if pad_x < 0:
154     raise ValueError("setFractionalPadding: Arguments must be non-negative")
155     if pad_x > 10:
156     raise ValueError("setFractionalPadding: Argument too large")
157     if pad_y is not None:
158     if pad_y < 0:
159     raise ValueError("setFractionalPadding: Arguments must be non-negative")
160     if pad_y > 10:
161     raise ValueError("setFractionalPadding: Argument too large")
162     self._padding = [pad_x,pad_y], 'f'
163 caltinay 4213
164 gross 4376 def setPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
165 caltinay 4108 """
166     Sets the amount of padding around the dataset in absolute length units.
167 caltinay 4060
168 caltinay 4108 The final domain size will be the length in x (in y) of the dataset
169     plus twice the value of `pad_x` (`pad_y`). The arguments must be
170     non-negative.
171    
172     :param pad_x: Padding per side in x direction (default: no padding)
173 gross 4373 :type pad_x: ``float`` in units of length (meter)
174 caltinay 4108 :param pad_y: Padding per side in y direction (default: no padding)
175 gross 4373 :type pad_y: ``float`` in units of length (meter)
176 caltinay 4150 :note: `pad_y` is ignored for 2-dimensional domains.
177 gross 4373 :note: this function can only be used if the reference system is Cartesian
178 caltinay 4108 """
179 gross 4373 if not self.getReferenceSystem().isCartesian():
180 caltinay 4394 raise RuntimeError("setPadding can be called for the Cartesian reference system only.")
181 caltinay 4150 if self.__domain is not None:
182     raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
183 caltinay 4108 if pad_x is not None:
184     if pad_x < 0:
185     raise ValueError("setPadding: Arguments must be non-negative")
186     if pad_y is not None:
187     if pad_y < 0:
188     raise ValueError("setPadding: Arguments must be non-negative")
189     self._padding = [pad_x,pad_y], 'l'
190 gross 4373
191     def setGeoPadding(self, pad_lat=None, pad_lon=None):
192     """
193     Sets the amount of padding around the dataset in longitude and latitude.
194 caltinay 4108
195 caltinay 4435 The final domain size will be the extent in the latitudinal (in
196     longitudinal) direction of the dataset plus twice the value of
197     `pad_lat` (`pad_lon`). The arguments must be non-negative.
198 gross 4373
199 caltinay 4435 :param pad_lat: Padding per side in latitudinal direction (default: 0)
200 gross 4373 :type pad_lat: ``float`` in units of degree
201 caltinay 4435 :param pad_lon: Padding per side in longitudinal direction (default: 0)
202 gross 4373 :type pad_lon: ``float`` in units of degree
203     :note: `pad_lon` is ignored for 2-dimensional domains.
204     :note: this function can only be used if the reference system is not Cartesian
205 caltinay 4108 """
206 gross 4373 if self.getReferenceSystem().isCartesian():
207 caltinay 4394 raise RuntimeError("setGeoPadding can be called for non-Cartesian reference systems only.")
208 gross 4373 if self.__domain is not None:
209     raise RuntimeError("Invalid call to setPadding(). Domain is already built.")
210     if pad_lat is not None:
211     if pad_lat < 0:
212     raise ValueError("setPadding: Arguments must be non-negative")
213     if pad_lon is not None:
214     if pad_lon < 0:
215     raise ValueError("setPadding: Arguments must be non-negative")
216     self._padding = [pad_lat,pad_lon], 'd'
217    
218     def setElementPadding(self, pad_x=None, pad_y=None, pad_lat=None, pad_lon=None):
219     """
220 caltinay 4145 Sets the amount of padding around the dataset in number of elements
221     (cells).
222 caltinay 4108
223     When the domain is constructed `pad_x` (`pad_y`) elements are added
224     on each side of the x- (y-) dimension. The arguments must be
225     non-negative.
226    
227     :param pad_x: Padding per side in x direction (default: no padding)
228     :type pad_x: ``int``
229     :param pad_y: Padding per side in y direction (default: no padding)
230     :type pad_y: ``int``
231     :note: `pad_y` is ignored for 2-dimensional datasets.
232     """
233 gross 4373 if not pad_lat == None:
234 caltinay 4394 if not pad_x == None:
235     raise ValueError("Either pad_lat or pad_x can be set.")
236     else:
237     pad_x = pad_lat
238 gross 4373 if not pad_lon == None:
239 caltinay 4394 if not pad_y == None:
240     raise ValueError("Either pad_lon or pad_y can be set.")
241     else:
242     pad_y = pad_lan
243    
244 caltinay 4150 if self.__domain is not None:
245     raise RuntimeError("Invalid call to setElementPadding(). Domain is already built.")
246 caltinay 4108 if pad_x is not None:
247     if type(pad_x) is not int:
248     raise TypeError("setElementPadding expects integer arguments")
249     if pad_x < 0:
250     raise ValueError("setElementPadding: Arguments must be non-negative")
251     if pad_y is not None:
252     if type(pad_y) is not int:
253     raise TypeError("setElementPadding expects integer arguments")
254     if pad_y < 0:
255     raise ValueError("setElementPadding: Arguments must be non-negative")
256     self._padding = [pad_x,pad_y], 'e'
257 caltinay 4213
258 caltinay 4060 def getGravitySurveys(self):
259 caltinay 4145 """
260     Returns a list of gravity surveys, see `getSurveys` for details.
261     """
262     return self.getSurveys(DataSource.GRAVITY)
263 caltinay 4213
264 gross 4115 def getMagneticSurveys(self):
265 gross 4120 """
266 caltinay 4145 Returns a list of magnetic surveys, see `getSurveys` for details.
267 gross 4120 """
268 caltinay 4145 return self.getSurveys(DataSource.MAGNETIC)
269 caltinay 4213
270 caltinay 4145 def fixDensityBelow(self, depth=None):
271     """
272 caltinay 4435 Defines the depth below which the density anomaly is set to a given
273     value. If no value is given zero is assumed.
274 gross 4343
275 caltinay 4435 :param depth: depth below which the density is fixed. If not set, no
276     constraint at depth is applied.
277 gross 4343 :type depth: ``float``
278 caltinay 4145 """
279 gross 4120 self.__fix_density_below=depth
280 caltinay 4213
281 caltinay 4145 def fixSusceptibilityBelow(self, depth=None):
282 gross 4120 """
283 caltinay 4435 Defines the depth below which the susceptibility anomaly is set to a
284     given value. If no value is given zero is assumed.
285 gross 4343
286 caltinay 4435 :param depth: depth below which the susceptibility is fixed. If not
287     set, no constraint at depth is applied.
288 gross 4343 :type depth: ``float``
289 gross 4120 """
290     self.__fix_susceptibility_below=depth
291    
292 gross 4115 def getSurveys(self, datatype):
293 caltinay 4060 """
294 caltinay 4145 Returns a list of `Data` objects for all surveys of type `datatype`
295     available to this domain builder.
296 caltinay 4060
297 caltinay 4150 :return: List of surveys which are tuples (anomaly,error).
298 caltinay 4060 :rtype: ``list``
299     """
300 gross 4124 surveys=[]
301 caltinay 4145 for src in self.__sources:
302     if src.getDataType()==datatype:
303     surveys.append(src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing))
304 gross 4124 return surveys
305 caltinay 4060
306 gross 4115 def setBackgroundMagneticFluxDensity(self, B):
307 caltinay 4060 """
308 gross 4433 Sets the background magnetic flux density B=(B_East, B_North, B_Vertical)
309 caltinay 4060 """
310 gross 4106 self.__background_magnetic_field=B
311 caltinay 4108
312 gross 4115 def getBackgroundMagneticFluxDensity(self):
313 gross 4106 """
314 caltinay 4145 Returns the background magnetic flux density.
315 caltinay 4108 """
316 gross 4106 B = self.__background_magnetic_field
317 caltinay 4462 if B is None:
318     raise ValueError("No background magnetic flux density set!")
319    
320 gross 4433 if self.__dim < 3 :
321 caltinay 4287 return np.array([B[0], B[2]])
322 caltinay 4060 else:
323 caltinay 4287 return np.array(B)
324 caltinay 4060
325     def getSetDensityMask(self):
326 caltinay 4150 """
327     Returns the density mask data object which is non-zero for cells
328     whose density value is fixed, zero otherwise.
329     """
330 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
331 gross 4120 m = whereNonNegative(z)
332     if self.__fix_density_below:
333 caltinay 4435 m += whereNonPositive(z+self.__v_scale*self.__fix_density_below)
334 gross 4120 return m
335 caltinay 4145
336 caltinay 4060 def getSetSusceptibilityMask(self):
337 caltinay 4150 """
338     Returns the susceptibility mask data object which is non-zero for
339     cells whose susceptibility value is fixed, zero otherwise.
340     """
341 caltinay 4145 z=self.getDomain().getX()[self.__dim-1]
342 gross 4120 m = whereNonNegative(z)
343     if self.__fix_susceptibility_below:
344 caltinay 4435 m += whereNonPositive(z+self.__v_scale*self.__fix_susceptibility_below)
345 gross 4120 return m
346 caltinay 4060
347     def getDomain(self):
348     """
349     Returns a domain that spans the data area plus padding.
350    
351 caltinay 4145 The domain is created the first time this method is called,
352     subsequent calls return the same domain so anything that affects
353     the domain (such as padding) needs to be set beforehand.
354    
355     :return: The escript domain for this data source
356 caltinay 4060 :rtype: `esys.escript.Domain`
357     """
358 caltinay 4145 if self.__domain is None:
359     self.__domain=self.__createDomain()
360     return self.__domain
361 caltinay 4060
362     def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
363     """
364     This method sets the target domain parameters for the vertical
365     dimension.
366    
367 caltinay 4145 :param depth: Depth of the domain (in meters)
368 caltinay 4060 :type depth: ``float``
369 caltinay 4145 :param air_layer: Depth of the layer above sea level (in meters)
370 caltinay 4060 :type air_layer: ``float``
371     :param num_cells: Number of domain elements for the entire vertical
372     dimension
373     :type num_cells: ``int``
374     """
375 caltinay 4150 if self.__domain is not None:
376     raise RuntimeError("Invalid call to setVerticalExtents(). Domain is already built.")
377 caltinay 4060 self._v_depth=depth
378     self._v_air_layer=air_layer
379     self._v_num_cells=num_cells
380    
381     def __getTotalExtentsWithPadding(self):
382     """
383 caltinay 4362 Helper method that computes origin and number of data elements
384 caltinay 4060 after adding padding to the bounding box of all available survey data.
385     """
386     X0, NX, DX = self.__getTotalExtents()
387 caltinay 4362 DATA_DIM=len(X0)
388 caltinay 4060 frac=[]
389     # padding is applied to each side so multiply by 2 to get the total
390     # amount of padding per dimension
391 caltinay 4108 pad, pt = self._padding
392 caltinay 4362 for i in range(DATA_DIM):
393 azadeh 4127 if pad[i] is None:
394 caltinay 4145 frac.append(0.)
395     continue
396 gross 4373 if pt == 'f' : # fraction of side length
397 caltinay 4108 frac.append(2.*pad[i])
398     elif pt == 'e': # number of elements
399     frac.append(2.*pad[i]/float(NX[i]))
400     else: # absolute length
401     f=pad[i]/DX[i]
402     frac.append(2.*f/float(NX[i]))
403 caltinay 4060
404     # calculate new number of elements
405 caltinay 4362 NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DATA_DIM)]
406     NXdiff=[NX_padded[i]-NX[i] for i in range(DATA_DIM)]
407     X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DATA_DIM)]
408 caltinay 4060 return X0_padded, NX_padded, DX
409    
410     def __getTotalExtents(self):
411     """
412     Helper method that computes the origin, number of elements and
413     minimal element spacing taking into account all available survey data.
414     """
415 caltinay 4145 if len(self.__sources)==0:
416 caltinay 4060 raise ValueError("No data")
417 caltinay 4145 X0, NE, DX = self.__sources[0].getDataExtents()
418 caltinay 4613 # do not mess with the values if only one source used
419     if len(self.__sources)>1:
420     XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
421 caltinay 4060
422 caltinay 4613 for src in self.__sources[1:]:
423     d_x0, d_ne, d_dx = src.getDataExtents()
424     for i in range(len(d_x0)):
425     X0[i]=min(X0[i], d_x0[i])
426     for i in range(len(d_dx)):
427     DX[i]=min(DX[i], d_dx[i])
428     for i in range(len(d_ne)):
429     XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
430     # FIXME: should this be rounded up instead?
431     NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
432 caltinay 4060 return X0, NE, DX
433    
434     def __createDomain(self):
435     """
436     Creates and returns an escript domain that spans the entire area of
437     available data plus a padding zone. This method is called only once
438     the first time `getDomain()` is invoked.
439    
440     :return: The escript domain
441     :rtype: `esys.escript.Domain`
442     """
443     X0, NX, DX = self.__getTotalExtentsWithPadding()
444    
445     # number of domain elements
446     NE = NX + [self._v_num_cells]
447    
448     # origin of domain
449 caltinay 4435 origin = X0 + [-self._v_depth*self.__v_scale]
450    
451 caltinay 4394 if self.getReferenceSystem().isCartesian():
452     # rounding will give us about meter-accuracy with UTM coordinates
453     self._dom_origin = [np.floor(oi) for oi in origin]
454     else:
455     # this should give us about meter-accuracy with lat/lon coords
456     self._dom_origin = [1e-5*np.floor(oi*1e5) for oi in origin]
457 caltinay 4435
458 caltinay 4060 # cell size / point spacing
459 caltinay 4435 spacing = DX + [self.__v_scale*np.floor((self._v_depth+self._v_air_layer)/self._v_num_cells)]
460 gross 4373 #self._spacing = [float(np.floor(si)) for si in spacing]
461     self._spacing = spacing
462 gross 4433
463 caltinay 4435 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self.__dim)]
464 caltinay 4613
465 caltinay 4145 if self.__dim==3:
466 caltinay 4060 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
467     else:
468 caltinay 4434 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
469 caltinay 4108
470 caltinay 4060 # ripley may internally adjust NE and length, so recompute
471 caltinay 4145 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self.__dim)]
472 caltinay 4435 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self.__dim)]
473 caltinay 4060
474     self.logger.debug("Domain size: "+str(self._dom_NE))
475     self.logger.debug(" length: "+str(self._dom_len))
476     self.logger.debug(" origin: "+str(self._dom_origin))
477     return dom
478    

  ViewVC Help
Powered by ViewVC 1.1.26