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

Contents of /trunk/downunder/py_src/datasources.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4362 - (show annotations)
Tue Apr 16 04:37:13 2013 UTC (6 years, 4 months ago) by caltinay
File MIME type: text/x-python
File size: 45315 byte(s)
Implemented NumpyData downunder source (w/ tests).
Added check of data dimensionality to domainbuilder since data_dim must equal
domain_dim-1 at the moment.

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2013 by University of Queensland
5 # 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 """Data readers/providers for inversions"""
17
18 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
19 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__ = ['DataSource', 'ErMapperData', 'NumpyData',\
26 'SyntheticDataBase', 'SyntheticFeatureData', 'SyntheticData',
27 'SmoothAnomaly']
28
29 import logging
30 import numpy as np
31 import tempfile
32 from esys.escript import ReducedFunction
33 from esys.escript import unitsSI as U
34 from esys.escript.linearPDEs import LinearSinglePDE
35 from esys.escript.util import *
36 from esys.ripley import ripleycpp
37
38 try:
39 from scipy.io.netcdf import netcdf_file
40 __all__ += ['NetCdfData']
41 except:
42 pass
43
44 def LatLonToUTM(lon, lat, wkt_string=None):
45 """
46 Converts one or more longitude,latitude pairs to the corresponding x,y
47 coordinates in the Universal Transverse Mercator projection.
48
49 :note: The ``pyproj`` module is required for this method. If it is not
50 found an exception is raised.
51 :note: If `wkt_string` is not given or invalid or the ``gdal`` module is
52 not available to convert the string, then the input values are
53 assumed to be using the Clarke 1866 ellipsoid.
54
55 :param lon: longitude value(s)
56 :type lon: ``float``, ``list``, ``tuple``, or ``numpy.array``
57 :param lat: latitude value(s)
58 :type lat: ``float``, ``list``, ``tuple``, or ``numpy.array``
59 :param wkt_string: Well-known text (WKT) string describing the coordinate
60 system used. The ``gdal`` module is used to convert
61 the string to the corresponding Proj4 string.
62 :type wkt_string: ``str``
63 :rtype: ``tuple``
64 """
65
66 try:
67 import pyproj
68 except:
69 raise ImportError("In order to perform coordinate transformations on "
70 "the data you are using the 'pyproj' Python module is required but "
71 "was not found. Please install the module and try again.")
72
73 # determine UTM zone from the input data
74 zone=int(np.median((np.floor((np.array(lon) + 180)/6) + 1) % 60))
75 try:
76 import osgeo.osr
77 srs = osgeo.osr.SpatialReference()
78 srs.ImportFromWkt(wkt_string)
79 p_src = pyproj.Proj(srs.ExportToProj4())
80 except:
81 p_src = pyproj.Proj('+proj=longlat +ellps=clrk66 +no_defs')
82
83 # check for hemisphere
84 if np.median(np.array(lat)) < 0.:
85 south='+south '
86 else:
87 south=''
88 p_dest = pyproj.Proj('+proj=utm +zone=%d %s+units=m'%(zone,south))
89 x,y=pyproj.transform(p_src, p_dest, lon, lat)
90 return x,y,zone
91
92 class DataSource(object):
93 """
94 A class that provides survey data for the inversion process.
95 This is an abstract base class that implements common functionality.
96 Methods to be overwritten by subclasses are marked as such.
97 This class assumes 2D data which is mapped to a slice of a 3D domain.
98 For other setups override the methods as required.
99 """
100
101 GRAVITY, MAGNETIC = list(range(2))
102
103 def __init__(self):
104 """
105 Constructor. Sets some defaults and initializes logger.
106 """
107 self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
108 self.__subsampling_factor=1
109
110 def getDataExtents(self):
111 """
112 returns a tuple of tuples ``( (x0, y0), (nx, ny), (dx, dy) )``, where
113
114 - ``x0``, ``y0`` = coordinates of data origin
115 - ``nx``, ``ny`` = number of data points in x and y
116 - ``dx``, ``dy`` = spacing of data points in x and y
117
118 This method must be implemented in subclasses.
119 """
120 raise NotImplementedError
121
122 def getDataType(self):
123 """
124 Returns the type of survey data managed by this source.
125 Subclasses must return `GRAVITY` or `MAGNETIC` as appropriate.
126 """
127 raise NotImplementedError
128
129 def getSurveyData(self, domain, origin, NE, spacing):
130 """
131 This method is called by the `DomainBuilder` to retrieve the survey
132 data as `Data` objects on the given domain.
133
134 Subclasses should return one or more `Data` objects with survey data
135 interpolated on the given `escript` domain. The exact return type
136 depends on the type of data.
137
138 :param domain: the escript domain to use
139 :type domain: `esys.escript.Domain`
140 :param origin: the origin coordinates of the domain
141 :type origin: ``tuple`` or ``list``
142 :param NE: the number of domain elements in each dimension
143 :type NE: ``tuple`` or ``list``
144 :param spacing: the cell sizes (node spacing) in the domain
145 :type spacing: ``tuple`` or ``list``
146 """
147 raise NotImplementedError
148
149 def getUtmZone(self):
150 """
151 All data source coordinates are converted to UTM (Universal Transverse
152 Mercator) in order to have useful domain extents. Subclasses should
153 implement this method and return the UTM zone number of the projected
154 coordinates.
155 """
156 raise NotImplementedError
157
158 def setSubsamplingFactor(self, f):
159 """
160 Sets the data subsampling factor (default=1).
161
162 The factor is applied in all dimensions. For example a 2D dataset
163 with 300 x 150 data points will be reduced to 150 x 75 when a
164 subsampling factor of 2 is used.
165 This becomes important when adding data of varying resolution to
166 a `DomainBuilder`.
167 """
168 self.__subsampling_factor=f
169
170 def getSubsamplingFactor(self):
171 """
172 Returns the subsampling factor that was set via `setSubsamplingFactor`
173 (see there).
174 """
175 return self.__subsampling_factor
176
177
178 ##############################################################################
179 class ErMapperData(DataSource):
180 """
181 Data Source for ER Mapper raster data.
182 Note that this class only accepts a very specific type of ER Mapper data
183 input and will raise an exception if other data is found.
184 """
185 def __init__(self, data_type, headerfile, datafile=None, altitude=0.,
186 error=None, scale_factor=None, null_value=None):
187 """
188 :param data_type: type of data, must be `GRAVITY` or `MAGNETIC`
189 :type data_type: ``int``
190 :param headerfile: ER Mapper header file (usually ends in .ers)
191 :type headerfile: ``str``
192 :param datafile: ER Mapper binary data file name. If not supplied the
193 name of the header file without '.ers' is assumed
194 :type datafile: ``str``
195 :param altitude: altitude of measurements above ground in meters
196 :type altitude: ``float``
197 :param error: constant value to use for the data uncertainties.
198 If a value is supplied, it is scaled by the same factor
199 as the measurements. If not provided the error is
200 assumed to be 2 units for all measurements (i.e. 0.2
201 mGal and 2 nT for gravity and magnetic, respectively)
202 :type error: ``float``
203 :param scale_factor: the measurements and error values are scaled by
204 this factor. By default, gravity data is assumed
205 to be given in 1e-6 m/s^2 (0.1 mGal), while
206 magnetic data is assumed to be in 1e-9 T (1 nT).
207 :type scale_factor: ``float``
208 :param null_value: value that is used in the file to mark undefined
209 areas. This information is usually included in the
210 file.
211 :type null_value: ``float``
212 """
213 super(ErMapperData, self).__init__()
214 self.__headerfile=headerfile
215 if datafile is None:
216 self.__datafile=headerfile[:-4]
217 else:
218 self.__datafile=datafile
219 self.__altitude=altitude
220 self.__data_type=data_type
221 self.__utm_zone = None
222 self.__scale_factor = scale_factor
223 self.__null_value = null_value
224 self.__error_value = error
225 self.__readHeader()
226
227 def __readHeader(self):
228 self.logger.debug("Checking Data Source: %s (header: %s)"%(self.__datafile, self.__headerfile))
229 metadata=open(self.__headerfile, 'r').readlines()
230 start=-1
231 for i in range(len(metadata)):
232 if metadata[i].strip() == 'DatasetHeader Begin':
233 start=i+1
234 if start==-1:
235 raise RuntimeError('Invalid ER Mapper header file ("DatasetHeader" not found)')
236
237 # parse header file filling dictionary of found values
238 md_dict={}
239 section=[]
240 for i in range(start, len(metadata)):
241 line=metadata[i].strip()
242 if line[-6:].strip() == 'Begin':
243 section.append(line[:-6].strip())
244 elif line[-4:].strip() == 'End':
245 if len(section)>0:
246 section.pop()
247 else:
248 vals=line.split('=')
249 if len(vals)==2:
250 key = vals[0].strip()
251 value = vals[1].strip()
252 fullkey='.'.join(section+[key])
253 md_dict[fullkey]=value
254
255 # check that that the data format/type is supported
256 try:
257 if md_dict['ByteOrder'] != 'LSBFirst':
258 raise RuntimeError('Unsupported byte order '+md_dict['ByteOrder'])
259 except KeyError:
260 self.logger.warn("Byte order not specified. Assuming LSB first.")
261
262 try:
263 if md_dict['DataType'] != 'Raster':
264 raise RuntimeError('Unsupported data type '+md_dict['DataType'])
265 except KeyError:
266 self.logger.warn("Data type not specified. Assuming raster data.")
267
268 try:
269 if md_dict['RasterInfo.CellType'] != 'IEEE4ByteReal':
270 raise RuntimeError('Unsupported data type '+md_dict['RasterInfo.CellType'])
271 except KeyError:
272 self.logger.warn("Cell type not specified. Assuming IEEE4ByteReal.")
273
274 try:
275 fileOffset = int(md_dict['HeaderOffset'])
276 except:
277 fileOffset = 0
278 if fileOffset > 0:
279 raise RuntimeError("ER Mapper data with header offset >0 not supported.")
280
281 # now extract required information
282 try:
283 NX = int(md_dict['RasterInfo.NrOfCellsPerLine'])
284 NY = int(md_dict['RasterInfo.NrOfLines'])
285 except:
286 raise RuntimeError("Could not determine extents of data")
287
288 ### mask/null value
289 # note that NaN is always filtered out in ripley
290 if self.__null_value is None:
291 try:
292 self.__null_value = float(md_dict['RasterInfo.NullCellValue'])
293 except:
294 self.logger.debug("Could not determine null value, using default.")
295 self.__null_value = 99999
296 elif not isinstance(self.__null_value,float) and not isinstance(self.__null_value,int):
297 raise TypeError("Invalid type of null_value parameter")
298
299 try:
300 spacingX = float(md_dict['RasterInfo.CellInfo.Xdimension'])
301 spacingY = float(md_dict['RasterInfo.CellInfo.Ydimension'])
302 except:
303 raise RuntimeError("Could not determine cell dimensions")
304
305 try:
306 if md_dict['CoordinateSpace.CoordinateType']=='EN':
307 originX = float(md_dict['RasterInfo.RegistrationCoord.Eastings'])
308 originY = float(md_dict['RasterInfo.RegistrationCoord.Northings'])
309 elif md_dict['CoordinateSpace.CoordinateType']=='LL':
310 originX = float(md_dict['RasterInfo.RegistrationCoord.Longitude'])
311 originY = float(md_dict['RasterInfo.RegistrationCoord.Latitude'])
312 else:
313 raise RuntimeError("Unknown CoordinateType")
314 except:
315 self.logger.warn("Could not determine coordinate origin. Setting to (0.0, 0.0)")
316 originX,originY = 0.0, 0.0
317
318 # data sets have origin in top-left corner so y runs top-down and
319 # we need to flip accordingly
320 originY-=NY*spacingY
321
322 if 'GEODETIC' in md_dict['CoordinateSpace.Projection']:
323 # it appears we have lat/lon coordinates so need to convert
324 # origin and spacing. Try using gdal to get the wkt if available:
325 try:
326 from osgeo import gdal
327 ds=gdal.Open(self.__headerfile)
328 wkt=ds.GetProjection()
329 except:
330 wkt='GEOGCS["GEOCENTRIC DATUM of AUSTRALIA",DATUM["GDA94",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'
331 self.logger.warn('GDAL not available or file read error, assuming GDA94 data')
332 originX_UTM,originY_UTM,zone=LatLonToUTM(originX, originY, wkt)
333 op1X,op1Y,_=LatLonToUTM(originX+spacingX, originY+spacingY, wkt)
334 self.__utm_zone = zone
335 # we are rounding to avoid interpolation issues
336 spacingX=np.round(op1X-originX_UTM)
337 spacingY=np.round(op1Y-originY_UTM)
338 originX=np.round(originX_UTM)
339 originY=np.round(originY_UTM)
340
341 self.__dataorigin=[originX, originY]
342 self.__delta = [spacingX, spacingY]
343 self.__nPts = [NX, NY]
344 self.__origin = [originX, originY]
345 ### scale factor
346 if self.__scale_factor is None:
347 if self.__data_type == self.GRAVITY:
348 self.logger.info("Assuming gravity data scale is 1e-6 m/s^2.")
349 self.__scale_factor = 1e-6
350 else:
351 self.logger.info("Assuming magnetic data units are 'nT'.")
352 self.__scale_factor = 1e-9
353
354 ### error value
355 if self.__error_value is None:
356 self.__error_value = 2.
357 elif not isinstance(self.__error_value,float) and not isinstance(self.__error_value,int):
358 raise TypeError("Invalid type of error parameter")
359
360
361 def getDataExtents(self):
362 """
363 returns ( (x0, y0), (nx, ny), (dx, dy) )
364 """
365 return (list(self.__origin), list(self.__nPts), list(self.__delta))
366
367 def getDataType(self):
368 return self.__data_type
369
370 def getSurveyData(self, domain, origin, NE, spacing):
371 FS = ReducedFunction(domain)
372 nValues=self.__nPts
373 # determine base location of this dataset within the domain
374 first=[int((self.__origin[i]-origin[i])/spacing[i]) for i in range(len(self.__nPts))]
375 # determine the resolution difference between domain and data.
376 # If domain has twice the resolution we can double up the data etc.
377 multiplier=[int(round(self.__delta[i]/spacing[i])) for i in range(len(self.__nPts))]
378 if domain.getDim()==3:
379 first.append(int((self.__altitude-origin[2])/spacing[2]))
380 multiplier=multiplier+[1]
381 nValues=nValues+[1]
382
383 data = ripleycpp._readBinaryGrid(self.__datafile, FS, first, nValues,
384 multiplier, (), self.__null_value)
385 sigma = self.__error_value * whereNonZero(data-self.__null_value)
386
387 data = data * self.__scale_factor
388 sigma = sigma * self.__scale_factor
389 return data, sigma
390
391 def getUtmZone(self):
392 return self.__utm_zone
393
394
395 ##############################################################################
396 class NetCdfData(DataSource):
397 """
398 Data Source for gridded netCDF data that use CF/COARDS conventions.
399 """
400 def __init__(self, data_type, filename, altitude=0., data_variable=None,
401 error=None, scale_factor=None, null_value=None):
402 """
403 :param filename: file name for survey data in netCDF format
404 :type filename: ``str``
405 :param data_type: type of data, must be `GRAVITY` or `MAGNETIC`
406 :type data_type: ``int``
407 :param altitude: altitude of measurements in meters
408 :type altitude: ``float``
409 :param data_variable: name of the netCDF variable that holds the data.
410 If not provided an attempt is made to determine
411 the variable and an exception thrown on failure.
412 :type data_variable: ``str``
413 :param error: either the name of the netCDF variable that holds the
414 uncertainties of the measurements or a constant value
415 to use for the uncertainties. If a constant value is
416 supplied, it is scaled by the same factor as the
417 measurements. If not provided the error is assumed to
418 be 2 units for all measurements (i.e. 0.2 mGal and 2 nT
419 for gravity and magnetic, respectively)
420 :type error: ``str`` or ``float``
421 :param scale_factor: the measurements and error values are scaled by
422 this factor. By default, gravity data is assumed
423 to be given in 1e-6 m/s^2 (0.1 mGal), while
424 magnetic data is assumed to be in 1e-9 T (1 nT).
425 :type scale_factor: ``float``
426 :param null_value: value that is used in the file to mark undefined
427 areas. This information is usually included in the
428 file.
429 :type null_value: ``float``
430 """
431 super(NetCdfData,self).__init__()
432 self.__filename=filename
433 if not data_type in [self.GRAVITY,self.MAGNETIC]:
434 raise ValueError("Invalid value for data_type parameter")
435 self.__data_type = data_type
436 self.__altitude = altitude
437 self.__data_name = data_variable
438 self.__scale_factor = scale_factor
439 self.__null_value = null_value
440 self.__readMetadata(error)
441
442 def __readMetadata(self, error):
443 self.logger.debug("Checking Data Source: %s"%self.__filename)
444 f=netcdf_file(self.__filename, 'r')
445
446 ### longitude- / X-dimension and variable
447 NX=0
448 for n in ['lon','longitude','x']:
449 if n in f.dimensions:
450 NX=f.dimensions[n]
451 lon_name=n
452 break
453 if NX==0:
454 raise RuntimeError("Could not determine extents of data")
455
456 # CF/COARDS states that coordinate variables have the same name as
457 # the dimensions
458 if not lon_name in f.variables:
459 raise RuntimeError("Could not determine longitude variable")
460 longitude=f.variables.pop(lon_name)
461
462 ### latitude- / Y-dimension and variable
463 NY=0
464 for n in ['lat','latitude','y']:
465 if n in f.dimensions:
466 NY=f.dimensions[n]
467 lat_name=n
468 break
469 if NY==0:
470 raise RuntimeError("Could not determine extents of data")
471 if not lat_name in f.variables:
472 raise RuntimeError("Could not determine latitude variable")
473 latitude=f.variables.pop(lat_name)
474
475 ### data variable
476 if self.__data_name is not None:
477 try:
478 dims = f.variables[self.__data_name].dimensions
479 if not ((lat_name in dims) and (lon_name in dims)):
480 raise ValueError("Invalid data variable name supplied")
481 except KeyError:
482 raise ValueError("Invalid data variable name supplied")
483 else:
484 for n in f.variables.keys():
485 dims=f.variables[n].dimensions
486 if (lat_name in dims) and (lon_name in dims):
487 self.__data_name=n
488 break
489 if self.__data_name is None:
490 raise RuntimeError("Could not determine data variable")
491
492 datavar = f.variables[self.__data_name]
493
494 ### error value/variable
495 self.__error_name = None
496 if isinstance(error,str):
497 try:
498 dims = f.variables[error].dimensions
499 if not ((lat_name in dims) and (lon_name in dims)):
500 raise ValueError("Invalid error variable name supplied")
501 except KeyError:
502 raise ValueError("Invalid error variable name supplied")
503 self.__error_name = error
504 elif isinstance(error,float) or isinstance(error,int):
505 self.__error_value = float(error)
506 elif error is None:
507 self.__error_value = 2.
508 else:
509 raise TypeError("Invalid type of error parameter")
510
511 ### mask/null value
512 # note that NaN is always filtered out in ripley
513 if self.__null_value is None:
514 if hasattr(datavar, 'missing_value'):
515 self.__null_value = float(datavar.missing_value)
516 elif hasattr(datavar, '_FillValue'):
517 self.__null_value = float(datavar._FillValue)
518 else:
519 self.logger.debug("Could not determine null value, using default.")
520 self.__null_value = 99999
521 elif not isinstance(self.__null_value,float) and not isinstance(self.__null_value,int):
522 raise TypeError("Invalid type of null_value parameter")
523
524 # try to determine units of data - this is disabled until we obtain a
525 # file with valid information
526 #if hasattr(f.variables[data_name], 'units'):
527 # units=f.variables[data_name].units
528
529 ### scale factor
530 if self.__scale_factor is None:
531 if self.__data_type == self.GRAVITY:
532 self.logger.info("Assuming gravity data scale is 1e-6 m/s^2.")
533 self.__scale_factor = 1e-6
534 else:
535 self.logger.info("Assuming magnetic data units are 'nT'.")
536 self.__scale_factor = 1e-9
537
538 # see if there is a WKT string to convert coordinates
539 try:
540 wkt_string=datavar.esri_pe_string
541 except:
542 wkt_string=None
543
544 # actual_range & geospatial_lon_min/max do not always contain correct
545 # values so we have to obtain the min/max in a less efficient way:
546 lon_range=longitude.data.min(),longitude.data.max()
547 lat_range=latitude.data.min(),latitude.data.max()
548 lon_range,lat_range,zone=LatLonToUTM(lon_range, lat_range, wkt_string)
549 self.__utm_zone = zone
550 lengths=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]
551 f.close()
552
553 self.__nPts=[NX, NY]
554 self.__origin=[lon_range[0],lat_range[0]]
555 # we are rounding to avoid interpolation issues
556 self.__delta=[np.round(lengths[i]/self.__nPts[i]) for i in range(2)]
557
558 def getDataExtents(self):
559 """
560 returns ( (x0, y0), (nx, ny), (dx, dy) )
561 """
562 return (list(self.__origin), list(self.__nPts), list(self.__delta))
563
564 def getDataType(self):
565 return self.__data_type
566
567 def getSurveyData(self, domain, origin, NE, spacing):
568 FS=ReducedFunction(domain)
569 nValues=self.__nPts
570 # determine base location of this dataset within the domain
571 first=[int((self.__origin[i]-origin[i])/spacing[i]) for i in range(len(self.__nPts))]
572 # determine the resolution difference between domain and data.
573 # If domain has twice the resolution we can double up the data etc.
574 multiplier=[int(round(self.__delta[i]/spacing[i])) for i in range(len(self.__nPts))]
575 if domain.getDim()==3:
576 first.append(int((self.__altitude-origin[2])/spacing[2]))
577 multiplier=multiplier+[1]
578 nValues=nValues+[1]
579
580 data = ripleycpp._readNcGrid(self.__filename, self.__data_name, FS,
581 first, nValues, multiplier, (),
582 self.__null_value)
583 if self.__error_name is not None:
584 sigma = ripleycpp._readNcGrid(self.__filename, self.__error_name,
585 FS, first, nValues, multiplier, (),0.)
586 else:
587 sigma = self.__error_value * whereNonZero(data-self.__null_value)
588
589 data = data * self.__scale_factor
590 sigma = sigma * self.__scale_factor
591 return data, sigma
592
593 def getUtmZone(self):
594 return self.__utm_zone
595
596
597 ##############################################################################
598 class SourceFeature(object):
599 """
600 A feature adds a density/susceptibility distribution to (parts of) a
601 domain of a synthetic data source, for example a layer of a specific
602 rock type or a simulated ore body.
603 """
604 def getValue(self):
605 """
606 Returns the value for the area covered by mask. It can be constant
607 or a `Data` object with spatial dependency.
608 """
609 raise NotImplementedError
610
611 def getMask(self, x):
612 """
613 Returns the mask of the area of interest for this feature. That is,
614 mask is non-zero where the value returned by `getValue()` should be
615 applied, zero elsewhere.
616 """
617 raise NotImplementedError
618
619 class SmoothAnomaly(SourceFeature):
620 """
621 A source feature in the form of a blob (roughly gaussian).
622 """
623 def __init__(self, lx, ly, lz, x, y, depth, v_inner=None, v_outer=None):
624 """
625 Intializes the smooth anomaly data.
626
627 :param lx: size of blob in x-dimension
628 :param ly: size of blob in y-dimension
629 :param lz: size of blob in z-dimension
630 :param x: location of blob in x-dimension
631 :param y: location of blob in y-dimension
632 :param depth: depth of blob
633 :param v_inner: value in the centre of the blob
634 :param v_outer: value in the periphery of the blob
635 """
636 self.x=x
637 self.y=y
638 self.lx=lx
639 self.ly=ly
640 self.lz=lz
641 self.depth=depth
642 self.v_inner=v_inner
643 self.v_outer=v_outer
644 self.value=None
645 self.mask=None
646
647 def getValue(self,x):
648 if self.value is None:
649 if self.v_outer is None or self.v_inner is None:
650 self.value=0
651 else:
652 DIM=x.getDomain().getDim()
653 alpha=-log(abs(self.v_outer/self.v_inner))*4
654 value=exp(-alpha*((x[0]-self.x)/self.lx)**2)
655 value=value*exp(-alpha*((x[DIM-1]+self.depth)/self.lz)**2)
656 self.value=maximum(abs(self.v_outer), abs(self.v_inner*value))
657 if self.v_inner<0: self.value=-self.value
658
659 return self.value
660
661 def getMask(self, x):
662 DIM=x.getDomain().getDim()
663 m=whereNonNegative(x[DIM-1]+self.depth+self.lz/2) * whereNonPositive(x[DIM-1]+self.depth-self.lz/2) \
664 *whereNonNegative(x[0]-(self.x-self.lx/2)) * whereNonPositive(x[0]-(self.x+self.lx/2))
665 if DIM>2:
666 m*=whereNonNegative(x[1]-(self.y-self.ly/2)) * whereNonPositive(x[1]-(self.y+self.ly/2))
667 self.mask = m
668 return m
669
670 ##############################################################################
671 class SyntheticDataBase(DataSource):
672 """
673 Base class to define reference data based on a given property distribution
674 (density or susceptibility). Data are collected from a square region of
675 vertical extent ``length`` on a grid with ``number_of_elements`` cells in
676 each direction.
677
678 The synthetic data are constructed by solving the appropriate forward
679 problem. Data can be sampled with an offset from the surface at z=0 or
680 using the entire subsurface region.
681 """
682 def __init__(self, data_type,
683 DIM=2,
684 number_of_elements=10,
685 length=1*U.km,
686 B_b=None,
687 data_offset=0,
688 full_knowledge=False,
689 spherical=False):
690 """
691 :param data_type: data type indicator
692 :type data_type: `DataSource.GRAVITY`, `DataSource.MAGNETIC`
693 :param DIM: number of spatial dimensions
694 :type DIM: ``int`` (2 or 3)
695 :param number_of_elements: lateral number of elements in the region
696 where data are collected
697 :type number_of_elements: ``int``
698 :param length: lateral extent of the region where data are collected
699 :type length: ``float``
700 :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude]. Only used for magnetic data.
701 :type B_b: ``list`` of ``Scalar``
702 :param data_offset: offset of the data collection region from the surface
703 :type data_offset: ``float``
704 :param full_knowledge: if ``True`` data are collected from the entire
705 subsurface region. This is mainly for testing.
706 :type full_knowledge: ``Bool``
707 :param spherical: if ``True`` spherical coordinates are used. This is
708 not supported yet and should be set to ``False``.
709 :type spherical: ``Bool``
710 """
711 super(SyntheticDataBase, self).__init__()
712 if not data_type in [self.GRAVITY, self.MAGNETIC]:
713 raise ValueError("Invalid value for data_type parameter")
714 self.DIM=DIM
715 self.number_of_elements=number_of_elements
716 self.length=length
717 self.__data_type = data_type
718
719 if spherical:
720 raise ValueError("Spherical coordinates are not supported yet")
721 self.__spherical = spherical
722 self.__full_knowledge= full_knowledge
723 self.__data_offset=data_offset
724 self.__B_b =None
725 # this is for Cartesian (FIXME ?)
726 if data_type == self.MAGNETIC:
727 if self.DIM < 3:
728 self.__B_b = np.array([-B_b[2], -B_b[0]])
729 else:
730 self.__B_b = ([-B_b[1], -B_b[2], -B_b[0]])
731 self.__origin = [0]*(DIM-1)
732 self.__delta = [float(length)/number_of_elements]*(DIM-1)
733 self.__nPts = [number_of_elements]*(DIM-1)
734 self._reference_data=None
735
736 def getUtmZone(self):
737 """
738 returns a dummy UTM zone since this class does not use real coordinate
739 values.
740 """
741 return 0
742
743 def getDataExtents(self):
744 """
745 returns the lateral data extend of the data set
746 """
747 return (list(self.__origin), list(self.__nPts), list(self.__delta))
748
749 def getDataType(self):
750 """
751 returns the data type
752 """
753 return self.__data_type
754
755 def getSurveyData(self, domain, origin, number_of_elements, spacing):
756 """
757 returns the survey data placed on a given domain.
758
759 :param domain: domain on which the data are to be placed
760 :type domain: ``Domain``
761 :param origin: origin of the domain
762 :type origin: ``list`` of ``float``
763 :param number_of_elements: number of elements (or cells) in each
764 spatial direction used to span the domain
765 :type number_of_elements: ``list`` of ``int``
766 :param spacing: cell size in each spatial direction
767 :type spacing: ``list`` of ``float``
768 :return: observed gravity field or magnetic flux density for each cell
769 in the domain and for each cell an indicator 1/0 if the data
770 are valid or not.
771 :rtype: pair of ``Scalar``
772 """
773 pde=LinearSinglePDE(domain)
774 DIM=domain.getDim()
775 x=domain.getX()
776 # set the reference data
777 k=self.getReferenceProperty(domain)
778 # calculate the corresponding potential
779 z=x[DIM-1]
780 m_psi_ref=whereZero(z-sup(z))
781 if self.getDataType()==DataSource.GRAVITY:
782 pde.setValue(A=kronecker(domain), Y=-4*np.pi*U.Gravitational_Constant*self._reference_data, q=m_psi_ref)
783 else:
784 pde.setValue(A=kronecker(domain), X=self._reference_data*self.__B_b, q=m_psi_ref)
785 pde.setSymmetryOn()
786 #pde.getSolverOptions().setTolerance(1e-13)
787 psi_ref=pde.getSolution()
788 del pde
789 if self.getDataType()==DataSource.GRAVITY:
790 data = -grad(psi_ref, ReducedFunction(domain))
791 else:
792 data = self._reference_data*self.__B_b-grad(psi_ref, ReducedFunction(domain))
793
794 x=ReducedFunction(domain).getX()
795 if self.__full_knowledge:
796 sigma = whereNegative(x[DIM-1])
797 else:
798 sigma=1.
799 # limit mask to non-padding in horizontal area
800 for i in range(DIM-1):
801 x_i=x[i]
802 sigma=sigma * wherePositive(x_i) * whereNegative(x_i-(sup(x_i)+inf(x_i)))
803 # limit mask to one cell thickness at z=0
804 z=x[DIM-1]
805 oo=int(self.__data_offset/spacing[DIM-1]+0.5)*spacing[DIM-1]
806 sigma = sigma * whereNonNegative(z-oo) * whereNonPositive(z-oo-spacing[DIM-1])
807 return data,sigma
808
809 def getReferenceProperty(self, domain=None):
810 """
811 Returns the reference `Data` object that was used to generate
812 the gravity/susceptibility anomaly data.
813
814 :return: the density or susceptibility anomaly used to create the
815 survey data
816 :note: it can be assumed that in the first call the ``domain``
817 argument is present so the actual anomaly data can be created.
818 In subsequent calls this may not be true.
819 :note: method needs to be overwritten
820 """
821 raise NotImplementedError()
822
823 class SyntheticFeatureData(SyntheticDataBase):
824 """
825 Uses a list of `SourceFeature` objects to define synthetic anomaly data.
826 """
827 def __init__(self, data_type,
828 features,
829 DIM=2,
830 number_of_elements=10,
831 length=1*U.km,
832 B_b=None,
833 data_offset=0,
834 full_knowledge=False,
835 spherical=False):
836 """
837 :param data_type: data type indicator
838 :type data_type: `DataSource.GRAVITY`, `DataSource.MAGNETIC`
839 :param features: list of features. It is recommended that the features
840 are located entirely below the surface.
841 :type features: ``list`` of `SourceFeature`
842 :param DIM: spatial dimensionality
843 :type DIM: ``int`` (2 or 3)
844 :param number_of_elements: lateral number of elements in the region
845 where data are collected
846 :type number_of_elements: ``int``
847 :param length: lateral extent of the region where data are collected
848 :type length: ``float``
849 :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude]. Only used for magnetic data.
850 :type B_b: ``list`` of ``Scalar``
851 :param data_offset: offset of the data collection region from the surface
852 :type data_offset: ``float``
853 :param full_knowledge: if ``True`` data are collected from the entire subsurface region. This is mainly for testing.
854 :type full_knowledge: ``Bool``
855 :param spherical: if ``True`` spherical coordinates are used.
856 :type spherical: ``Bool``
857 """
858 super(SyntheticFeatureData,self).__init__(
859 data_type=data_type, DIM=DIM,
860 number_of_elements=number_of_elements,
861 length=length, B_b=B_b,
862 data_offset=data_offset,
863 full_knowledge=full_knowledge,
864 spherical=spherical)
865 self._features = features
866
867
868 def getReferenceProperty(self, domain=None):
869 """
870 Returns the reference `Data` object that was used to generate
871 the gravity/susceptibility anomaly data.
872 """
873 if self._reference_data == None:
874 DIM=domain.getDim()
875 x=domain.getX()
876 k=0.
877 for f in self._features:
878 m=f.getMask(x)
879 k = k * (1-m) + f.getValue(x) * m
880 self._reference_data= k
881 return self._reference_data
882
883 class SyntheticData(SyntheticDataBase):
884 """
885 Defines synthetic gravity/magnetic data based on harmonic property anomaly
886
887 rho = amplitude * sin(n_depth * pi /depth * (z+depth_offset)) * sin(n_length * pi /length * (x - shift) )
888
889 for all x and z<=0. For z>0 rho = 0.
890 """
891 def __init__(self, data_type,
892 n_length=1,
893 n_depth=1,
894 depth_offset=0.,
895 depth=None,
896 amplitude=None,
897 DIM=2,
898 number_of_elements=10,
899 length=1*U.km,
900 B_b=None,
901 data_offset=0,
902 full_knowledge=False,
903 spherical=False,
904 s=0.):
905 """
906 :param data_type: data type indicator
907 :type data_type: `DataSource.GRAVITY`, `DataSource.MAGNETIC`
908 :param n_length: number of oscillations in the anomaly data within the
909 observation region
910 :type n_length: ``int``
911 :param n_depth: number of oscillations in the anomaly data below surface
912 :type n_depth: ``int``
913 :param depth_offset: vertical offset of the data
914 :type depth_offset: ``float``
915 :param depth: vertical extent in the anomaly data. If not present the
916 depth of the domain is used.
917 :type depth: ``float``
918 :param amplitude: data amplitude. Default value is 200 U.kg/U.m**3 for
919 gravity and 0.1 for magnetic data.
920 :param DIM: spatial dimensionality
921 :type DIM: ``int`` (2 or 3)
922 :param number_of_elements: lateral number of elements in the region
923 where data are collected
924 :type number_of_elements: ``int``
925 :param length: lateral extent of the region where data are collected
926 :type length: ``float``
927 :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude].
928 Only used for magnetic data.
929 :type B_b: ``list`` of ``Scalar``
930 :param data_offset: offset of the data collection region from the surface
931 :type data_offset: ``float``
932 :param full_knowledge: if ``True`` data are collected from the entire
933 subsurface region. This is mainly for testing.
934 :type full_knowledge: ``Bool``
935 :param spherical: if ``True`` spherical coordinates are used
936 :type spherical: ``Bool``
937 """
938 super(SyntheticData,self).__init__(
939 data_type=data_type, DIM=DIM,
940 number_of_elements=number_of_elements,
941 length=length, B_b=B_b,
942 data_offset=data_offset,
943 full_knowledge=full_knowledge,
944 spherical=spherical)
945 self.__n_length = n_length
946 self.__n_depth = n_depth
947 self.depth = depth
948 self.depth_offset=depth_offset
949 if amplitude == None:
950 if data_type == DataSource.GRAVITY:
951 amplitude = 200 *U.kg/U.m**3
952 else:
953 amplitude = 0.1
954 self.__amplitude = amplitude
955 self.__s=s
956
957 def getReferenceProperty(self, domain=None):
958 """
959 Returns the reference `Data` object that was used to generate
960 the gravity/susceptibility anomaly data.
961 """
962 if self._reference_data is None:
963 DIM=domain.getDim()
964 x=domain.getX()
965 # set the reference data
966 z=x[DIM-1]
967 dd=self.depth
968 if dd is None: dd=inf(z)
969 z2=(z+self.depth_offset)/(self.depth_offset-dd)
970 k=sin(self.__n_depth * np.pi * z2) * whereNonNegative(z2) * whereNonPositive(z2-1.) * self.__amplitude
971 for i in range(DIM-1):
972 x_i=x[i]
973 min_x=inf(x_i)
974 max_x=sup(x_i)
975 k*= sin(self.__n_length*np.pi*(x_i-min_x-self.__s)/(max_x-min_x))
976 self._reference_data= k
977 return self._reference_data
978
979 ##############################################################################
980 class NumpyData(DataSource):
981 """
982 """
983
984 def __init__(self, data_type, data, error=1., length=1.*U.km, null_value=-1.):
985 """
986 A data source that uses survey data from a ``numpy`` object or list
987 instead of a file.
988 The dimensionality is inferred from the shape of ``data`` (1- and
989 2-dimensional data is supported). The data origin is assumed to be
990 at the coordinate origin.
991
992 :param data_type: data type indicator
993 :type data_type: `DataSource.GRAVITY`, `DataSource.MAGNETIC`
994 :param data: the survey data array. Note that for a cartesian coordinate
995 system the shape of the data is considered to be
996 (nz,ny,nx).
997 :type data: ``numpy.array`` or ``list``
998 :param error: constant value to use for the data uncertainties or a
999 numpy object with uncertainties for every data point.
1000 :type error: ``float`` or ``list`` or ``ndarray``
1001 :param length: side length(s) of the data slice/volume. This can be
1002 a scalar to indicate constant length in all dimensions
1003 or an array/list of values in each coordinate dimension.
1004 :type length: ``float`` or ``list`` or ``ndarray``
1005 :param null_value: value that is used in the undefined regions of the
1006 survey data object.
1007 :type null_value: ``float``
1008 """
1009 super(NumpyData, self).__init__()
1010 if not data_type in [self.GRAVITY, self.MAGNETIC]:
1011 raise ValueError("Invalid value for data_type parameter")
1012 self.__data_type = data_type
1013 self.__data = np.asarray(data, dtype=np.float32)
1014 DIM = len(self.__data.shape)
1015 if DIM not in (1,2):
1016 raise ValueError("NumpyData requires 1- or 2-dimensional data")
1017 self.__error = np.asarray(error, dtype=np.float32)
1018 if len(self.__error.shape) > 0 and \
1019 self.__error.shape != self.__data.shape:
1020 raise ValueError("error argument must be scalar or match the shape of the data")
1021 if isinstance(length,float) or isinstance(length,int):
1022 length = [float(length)] * DIM
1023 else:
1024 length=np.asarray(length, dtype=float)
1025 if len(length.shape) != 1 or length.shape[0] != DIM:
1026 raise ValueError("length must be scalar or an array with %d values"%DIM)
1027 length=length.tolist()
1028
1029 self.__null_value = null_value
1030 self.__nPts = list(reversed(self.__data.shape))
1031 self.__delta = [length[i]/self.__nPts[i] for i in range(DIM)]
1032 self.__origin = [0.] * DIM # this could be an optional argument
1033
1034 def getDataExtents(self):
1035 return (self.__origin, self.__nPts, self.__delta)
1036
1037 def getDataType(self):
1038 return self.__data_type
1039
1040 def getSurveyData(self, domain, origin, NE, spacing):
1041 FS = ReducedFunction(domain)
1042 nValues = self.__nPts
1043 dataDIM = len(nValues)
1044 # determine base location of this dataset within the domain
1045 first=[int((self.__origin[i]-origin[i])/spacing[i]) for i in range(dataDIM)]
1046 # determine the resolution difference between domain and data.
1047 # If domain has twice the resolution we can double up the data etc.
1048 multiplier=[int(round(self.__delta[i]/spacing[i])) for i in range(dataDIM)]
1049 if domain.getDim() > dataDIM:
1050 first.append(int((-origin[-1])/spacing[-1]))
1051 multiplier=multiplier+[1]
1052 nValues=nValues+[1]
1053
1054 _handle, numpyfile = tempfile.mkstemp()
1055 os.close(_handle)
1056 self.__data.tofile(numpyfile)
1057
1058 data = ripleycpp._readBinaryGrid(numpyfile, FS, first, nValues,
1059 multiplier, (), self.__null_value)
1060 if len(self.__error.shape) > 0:
1061 self.__error.tofile(numpyfile)
1062 sigma = ripleycpp._readBinaryGrid(numpyfile, FS, first, nValues,
1063 multiplier, (), 0.)
1064 else:
1065 sigma = self.__error.item() * whereNonZero(data-self.__null_value)
1066
1067 os.unlink(numpyfile)
1068
1069 return data, sigma
1070
1071 def getUtmZone(self):
1072 """
1073 returns a dummy UTM zone since this class does not use real coordinate
1074 values.
1075 """
1076 return 0
1077

  ViewVC Help
Powered by ViewVC 1.1.26