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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.4120  
changed lines
  Added in v.4729

  ViewVC Help
Powered by ViewVC 1.1.26