/[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 4390 - (show annotations)
Mon May 6 01:59:34 2013 UTC (6 years, 5 months ago) by caltinay
File MIME type: text/x-python
File size: 46837 byte(s)
Fixed advertising of NetCdfData availability.

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

  ViewVC Help
Powered by ViewVC 1.1.26