/[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 4394 - (show annotations)
Tue May 7 04:56:59 2013 UTC (6 years, 5 months ago) by caltinay
File MIME type: text/x-python
File size: 46862 byte(s)
Fixed rounding with lat/lon coordinates and removed stray print statements.

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 op1X, op1Y= originX+spacingX, originY+spacingY
362 spacingX=np.round(op1X-originX,5)
363 spacingY=np.round(op1Y-originY,5)
364 originX=np.round(originX,5)
365 originY=np.round(originY,5)
366 self.__utm_zone = zone
367
368 self.__dataorigin=[originX, originY]
369 self.__delta = [spacingX, spacingY]
370 self.__nPts = [NX, NY]
371 self.__origin = [originX, originY]
372 ### scale factor
373 if self.__scale_factor is None:
374 if self.__data_type == self.GRAVITY:
375 self.logger.info("Assuming gravity data scale is 1e-6 m/s^2.")
376 self.__scale_factor = 1e-6
377 else:
378 self.logger.info("Assuming magnetic data units are 'nT'.")
379 self.__scale_factor = 1e-9
380
381 ### error value
382 if self.__error_value is None:
383 self.__error_value = 2.
384 elif not isinstance(self.__error_value,float) and not isinstance(self.__error_value,int):
385 raise TypeError("Invalid type of error parameter")
386
387
388 def getDataExtents(self):
389 """
390 returns ( (x0, y0), (nx, ny), (dx, dy) )
391 """
392 return (list(self.__origin), list(self.__nPts), list(self.__delta))
393
394 def getDataType(self):
395 return self.__data_type
396
397 def getSurveyData(self, domain, origin, NE, spacing):
398 FS = ReducedFunction(domain)
399 nValues=self.__nPts
400 # determine base location of this dataset within the domain
401 first=[int((self.__origin[i]-origin[i])/spacing[i]) for i in range(len(self.__nPts))]
402 # determine the resolution difference between domain and data.
403 # If domain has twice the resolution we can double up the data etc.
404 multiplier=[int(round(self.__delta[i]/spacing[i])) for i in range(len(self.__nPts))]
405 if domain.getDim()==3:
406 first.append(int((self.__altitude-origin[2])/spacing[2]))
407 multiplier=multiplier+[1]
408 nValues=nValues+[1]
409
410 data = ripleycpp._readBinaryGrid(self.__datafile, FS, first, nValues,
411 multiplier, (), self.__null_value)
412 sigma = self.__error_value * whereNonZero(data-self.__null_value)
413
414 data = data * self.__scale_factor
415 sigma = sigma * self.__scale_factor
416 return data, sigma
417
418 def getUtmZone(self):
419 return self.__utm_zone
420
421
422 ##############################################################################
423 class NetCdfData(DataSource):
424 """
425 Data Source for gridded netCDF data that use CF/COARDS conventions.
426 """
427 def __init__(self, data_type, filename, altitude=0., data_variable=None,
428 error=None, scale_factor=None, null_value=None, reference_system=None):
429 """
430 :param filename: file name for survey data in netCDF format
431 :type filename: ``str``
432 :param data_type: type of data, must be `GRAVITY` or `MAGNETIC`
433 :type data_type: ``int``
434 :param altitude: altitude of measurements in meters
435 :type altitude: ``float``
436 :param data_variable: name of the netCDF variable that holds the data.
437 If not provided an attempt is made to determine
438 the variable and an exception thrown on failure.
439 :type data_variable: ``str``
440 :param error: either the name of the netCDF variable that holds the
441 uncertainties of the measurements or a constant value
442 to use for the uncertainties. If a constant value is
443 supplied, it is scaled by the same factor as the
444 measurements. If not provided the error is assumed to
445 be 2 units for all measurements (i.e. 0.2 mGal and 2 nT
446 for gravity and magnetic, respectively)
447 :type error: ``str`` or ``float``
448 :param scale_factor: the measurements and error values are scaled by
449 this factor. By default, gravity data is assumed
450 to be given in 1e-6 m/s^2 (0.1 mGal), while
451 magnetic data is assumed to be in 1e-9 T (1 nT).
452 :type scale_factor: ``float``
453 :param null_value: value that is used in the file to mark undefined
454 areas. This information is usually included in the
455 file.
456 :type null_value: ``float``
457 :param reference_system: reference coordinate system to be used. For a Cartesian
458 reference the appropriate UTM transformation is applied.
459 By default the Cartesian coordinate system is used.
460 :type reference_system: `ReferenceSystem`
461 :note: consistence in the reference coordinate system and the reference coordinate system
462 used in the data source is not checked.
463 """
464 super(NetCdfData,self).__init__(reference_system)
465 self.__filename=filename
466 if not data_type in [self.GRAVITY,self.MAGNETIC]:
467 raise ValueError("Invalid value for data_type parameter")
468 self.__data_type = data_type
469 self.__altitude = altitude
470 self.__data_name = data_variable
471 self.__scale_factor = scale_factor
472 self.__null_value = null_value
473 self.__readMetadata(error)
474
475 def __readMetadata(self, error):
476 self.logger.debug("Checking Data Source: %s"%self.__filename)
477 f=netcdf_file(self.__filename, 'r')
478
479 ### longitude- / X-dimension and variable
480 NX=0
481 for n in ['lon','longitude','x']:
482 if n in f.dimensions:
483 NX=f.dimensions[n]
484 lon_name=n
485 break
486 if NX==0:
487 raise RuntimeError("Could not determine extents of data")
488
489 # CF/COARDS states that coordinate variables have the same name as
490 # the dimensions
491 if not lon_name in f.variables:
492 raise RuntimeError("Could not determine longitude variable")
493 longitude=f.variables.pop(lon_name)
494
495 ### latitude- / Y-dimension and variable
496 NY=0
497 for n in ['lat','latitude','y']:
498 if n in f.dimensions:
499 NY=f.dimensions[n]
500 lat_name=n
501 break
502 if NY==0:
503 raise RuntimeError("Could not determine extents of data")
504 if not lat_name in f.variables:
505 raise RuntimeError("Could not determine latitude variable")
506 latitude=f.variables.pop(lat_name)
507
508 ### data variable
509 if self.__data_name is not None:
510 try:
511 dims = f.variables[self.__data_name].dimensions
512 if not ((lat_name in dims) and (lon_name in dims)):
513 raise ValueError("Invalid data variable name supplied")
514 except KeyError:
515 raise ValueError("Invalid data variable name supplied")
516 else:
517 for n in f.variables.keys():
518 dims=f.variables[n].dimensions
519 if (lat_name in dims) and (lon_name in dims):
520 self.__data_name=n
521 break
522 if self.__data_name is None:
523 raise RuntimeError("Could not determine data variable")
524
525 datavar = f.variables[self.__data_name]
526
527 ### error value/variable
528 self.__error_name = None
529 if isinstance(error,str):
530 try:
531 dims = f.variables[error].dimensions
532 if not ((lat_name in dims) and (lon_name in dims)):
533 raise ValueError("Invalid error variable name supplied")
534 except KeyError:
535 raise ValueError("Invalid error variable name supplied")
536 self.__error_name = error
537 elif isinstance(error,float) or isinstance(error,int):
538 self.__error_value = float(error)
539 elif error is None:
540 self.__error_value = 2.
541 else:
542 raise TypeError("Invalid type of error parameter")
543
544 ### mask/null value
545 # note that NaN is always filtered out in ripley
546 if self.__null_value is None:
547 if hasattr(datavar, 'missing_value'):
548 self.__null_value = float(datavar.missing_value)
549 elif hasattr(datavar, '_FillValue'):
550 self.__null_value = float(datavar._FillValue)
551 else:
552 self.logger.debug("Could not determine null value, using default.")
553 self.__null_value = 99999
554 elif not isinstance(self.__null_value,float) and not isinstance(self.__null_value,int):
555 raise TypeError("Invalid type of null_value parameter")
556
557 # try to determine units of data - this is disabled until we obtain a
558 # file with valid information
559 #if hasattr(f.variables[data_name], 'units'):
560 # units=f.variables[data_name].units
561
562 ### scale factor
563 if self.__scale_factor is None:
564 if self.__data_type == self.GRAVITY:
565 self.logger.info("Assuming gravity data scale is 1e-6 m/s^2.")
566 self.__scale_factor = 1e-6
567 else:
568 self.logger.info("Assuming magnetic data units are 'nT'.")
569 self.__scale_factor = 1e-9
570
571 # see if there is a WKT string to convert coordinates
572 try:
573 wkt_string=datavar.esri_pe_string
574 except:
575 wkt_string=None
576
577 # actual_range & geospatial_lon_min/max do not always contain correct
578 # values so we have to obtain the min/max in a less efficient way:
579 lon_range=longitude.data.min(),longitude.data.max()
580 lat_range=latitude.data.min(),latitude.data.max()
581 if self.getReferenceSystem().isCartesian():
582 lon_range,lat_range,zone=LatLonToUTM(lon_range, lat_range, wkt_string)
583 else:
584 _,_,zone=LatLonToUTM(lon_range, lat_range, wkt_string)
585
586 self.__utm_zone = zone
587 lengths=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]
588 f.close()
589
590 self.__nPts=[NX, NY]
591 self.__origin=[lon_range[0],lat_range[0]]
592 # we are rounding to avoid interpolation issues
593 if self.getReferenceSystem().isCartesian():
594 # rounding will give us about meter-accuracy with UTM coordinates
595 r=0
596 else:
597 # this should give us about meter-accuracy with lat/lon coords
598 r=5
599 self.__delta=[np.round(lengths[i]/self.__nPts[i],r) for i in range(2)]
600
601 def getDataExtents(self):
602 """
603 returns ( (x0, y0), (nx, ny), (dx, dy) )
604 """
605 return (list(self.__origin), list(self.__nPts), list(self.__delta))
606
607 def getDataType(self):
608 return self.__data_type
609
610 def getSurveyData(self, domain, origin, NE, spacing):
611 FS=ReducedFunction(domain)
612 nValues=self.__nPts
613 # determine base location of this dataset within the domain
614 first=[int((self.__origin[i]-origin[i])/spacing[i]) for i in range(len(self.__nPts))]
615 # determine the resolution difference between domain and data.
616 # If domain has twice the resolution we can double up the data etc.
617 multiplier=[int(round(self.__delta[i]/spacing[i])) for i in range(len(self.__nPts))]
618 if domain.getDim()==3:
619 first.append(int((self.__altitude-origin[2])/spacing[2]))
620 multiplier=multiplier+[1]
621 nValues=nValues+[1]
622
623 data = ripleycpp._readNcGrid(self.__filename, self.__data_name, FS,
624 first, nValues, multiplier, (),
625 self.__null_value)
626 if self.__error_name is not None:
627 sigma = ripleycpp._readNcGrid(self.__filename, self.__error_name,
628 FS, first, nValues, multiplier, (),0.)
629 else:
630 sigma = self.__error_value * whereNonZero(data-self.__null_value)
631
632 data = data * self.__scale_factor
633 sigma = sigma * self.__scale_factor
634 return data, sigma
635
636 def getUtmZone(self):
637 return self.__utm_zone
638
639
640 ##############################################################################
641 class SourceFeature(object):
642 """
643 A feature adds a density/susceptibility distribution to (parts of) a
644 domain of a synthetic data source, for example a layer of a specific
645 rock type or a simulated ore body.
646 """
647 def getValue(self):
648 """
649 Returns the value for the area covered by mask. It can be constant
650 or a `Data` object with spatial dependency.
651 """
652 raise NotImplementedError
653
654 def getMask(self, x):
655 """
656 Returns the mask of the area of interest for this feature. That is,
657 mask is non-zero where the value returned by `getValue()` should be
658 applied, zero elsewhere.
659 """
660 raise NotImplementedError
661
662 class SmoothAnomaly(SourceFeature):
663 """
664 A source feature in the form of a blob (roughly gaussian).
665 """
666 def __init__(self, lx, ly, lz, x, y, depth, v_inner=None, v_outer=None):
667 """
668 Intializes the smooth anomaly data.
669
670 :param lx: size of blob in x-dimension
671 :param ly: size of blob in y-dimension
672 :param lz: size of blob in z-dimension
673 :param x: location of blob in x-dimension
674 :param y: location of blob in y-dimension
675 :param depth: depth of blob
676 :param v_inner: value in the centre of the blob
677 :param v_outer: value in the periphery of the blob
678 """
679 self.x=x
680 self.y=y
681 self.lx=lx
682 self.ly=ly
683 self.lz=lz
684 self.depth=depth
685 self.v_inner=v_inner
686 self.v_outer=v_outer
687 self.value=None
688 self.mask=None
689
690 def getValue(self,x):
691 if self.value is None:
692 if self.v_outer is None or self.v_inner is None:
693 self.value=0
694 else:
695 DIM=x.getDomain().getDim()
696 alpha=-log(abs(self.v_outer/self.v_inner))*4
697 value=exp(-alpha*((x[0]-self.x)/self.lx)**2)
698 value=value*exp(-alpha*((x[DIM-1]+self.depth)/self.lz)**2)
699 self.value=maximum(abs(self.v_outer), abs(self.v_inner*value))
700 if self.v_inner<0: self.value=-self.value
701
702 return self.value
703
704 def getMask(self, x):
705 DIM=x.getDomain().getDim()
706 m=whereNonNegative(x[DIM-1]+self.depth+self.lz/2) * whereNonPositive(x[DIM-1]+self.depth-self.lz/2) \
707 *whereNonNegative(x[0]-(self.x-self.lx/2)) * whereNonPositive(x[0]-(self.x+self.lx/2))
708 if DIM>2:
709 m*=whereNonNegative(x[1]-(self.y-self.ly/2)) * whereNonPositive(x[1]-(self.y+self.ly/2))
710 self.mask = m
711 return m
712
713 ##############################################################################
714 class SyntheticDataBase(DataSource):
715 """
716 Base class to define reference data based on a given property distribution
717 (density or susceptibility). Data are collected from a square region of
718 vertical extent ``length`` on a grid with ``number_of_elements`` cells in
719 each direction.
720
721 The synthetic data are constructed by solving the appropriate forward
722 problem. Data can be sampled with an offset from the surface at z=0 or
723 using the entire subsurface region.
724 """
725 def __init__(self, data_type,
726 DIM=2,
727 number_of_elements=10,
728 length=1*U.km,
729 B_b=None,
730 data_offset=0,
731 full_knowledge=False):
732 """
733 :param data_type: data type indicator
734 :type data_type: `DataSource.GRAVITY`, `DataSource.MAGNETIC`
735 :param DIM: number of spatial dimensions
736 :type DIM: ``int`` (2 or 3)
737 :param number_of_elements: lateral number of elements in the region
738 where data are collected
739 :type number_of_elements: ``int``
740 :param length: lateral extent of the region where data are collected
741 :type length: ``float``
742 :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude]. Only used for magnetic data.
743 :type B_b: ``list`` of ``Scalar``
744 :param data_offset: offset of the data collection region from the surface
745 :type data_offset: ``float``
746 :param full_knowledge: if ``True`` data are collected from the entire
747 subsurface region. This is mainly for testing.
748 :type full_knowledge: ``Bool``
749 """
750 super(SyntheticDataBase, self).__init__()
751 if not data_type in [self.GRAVITY, self.MAGNETIC]:
752 raise ValueError("Invalid value for data_type parameter")
753 self.DIM=DIM
754 self.number_of_elements=number_of_elements
755 self.length=length
756 self.__data_type = data_type
757 self.__full_knowledge= full_knowledge
758 self.__data_offset=data_offset
759 self.__B_b =None
760 # this is for Cartesian (FIXME ?)
761 if data_type == self.MAGNETIC:
762 if self.DIM < 3:
763 self.__B_b = np.array([B_b[0], B_b[2]])
764 else:
765 self.__B_b = ([B_b[0], B_b[1],B_b[2]])
766 self.__origin = [0]*(DIM-1)
767 self.__delta = [float(length)/number_of_elements]*(DIM-1)
768 self.__nPts = [number_of_elements]*(DIM-1)
769 self._reference_data=None
770
771 def getUtmZone(self):
772 """
773 returns a dummy UTM zone since this class does not use real coordinate
774 values.
775 """
776 return 0
777
778 def getDataExtents(self):
779 """
780 returns the lateral data extend of the data set
781 """
782 return (list(self.__origin), list(self.__nPts), list(self.__delta))
783
784 def getDataType(self):
785 """
786 returns the data type
787 """
788 return self.__data_type
789
790 def getSurveyData(self, domain, origin, number_of_elements, spacing):
791 """
792 returns the survey data placed on a given domain.
793
794 :param domain: domain on which the data are to be placed
795 :type domain: ``Domain``
796 :param origin: origin of the domain
797 :type origin: ``list`` of ``float``
798 :param number_of_elements: number of elements (or cells) in each
799 spatial direction used to span the domain
800 :type number_of_elements: ``list`` of ``int``
801 :param spacing: cell size in each spatial direction
802 :type spacing: ``list`` of ``float``
803 :return: observed gravity field or magnetic flux density for each cell
804 in the domain and for each cell an indicator 1/0 if the data
805 are valid or not.
806 :rtype: pair of ``Scalar``
807 """
808 pde=LinearSinglePDE(domain)
809 DIM=domain.getDim()
810 x=domain.getX()
811 # set the reference data
812 k=self.getReferenceProperty(domain)
813 # calculate the corresponding potential
814 z=x[DIM-1]
815 m_psi_ref=whereZero(z-sup(z))
816 if self.getDataType()==DataSource.GRAVITY:
817 pde.setValue(A=kronecker(domain), Y=-4*np.pi*U.Gravitational_Constant*self._reference_data, q=m_psi_ref)
818 else:
819 pde.setValue(A=kronecker(domain), X=self._reference_data*self.__B_b, q=m_psi_ref)
820 pde.setSymmetryOn()
821 #pde.getSolverOptions().setTolerance(1e-13)
822 psi_ref=pde.getSolution()
823 del pde
824 if self.getDataType()==DataSource.GRAVITY:
825 data = -grad(psi_ref, ReducedFunction(domain))
826 else:
827 data = self._reference_data*self.__B_b-grad(psi_ref, ReducedFunction(domain))
828
829 x=ReducedFunction(domain).getX()
830 if self.__full_knowledge:
831 sigma = whereNegative(x[DIM-1])
832 else:
833 sigma=1.
834 # limit mask to non-padding in horizontal area
835 for i in range(DIM-1):
836 x_i=x[i]
837 sigma=sigma * wherePositive(x_i) * whereNegative(x_i-(sup(x_i)+inf(x_i)))
838 # limit mask to one cell thickness at z=0
839 z=x[DIM-1]
840 oo=int(self.__data_offset/spacing[DIM-1]+0.5)*spacing[DIM-1]
841 sigma = sigma * whereNonNegative(z-oo) * whereNonPositive(z-oo-spacing[DIM-1])
842 return data,sigma
843
844 def getReferenceProperty(self, domain=None):
845 """
846 Returns the reference `Data` object that was used to generate
847 the gravity/susceptibility anomaly data.
848
849 :return: the density or susceptibility anomaly used to create the
850 survey data
851 :note: it can be assumed that in the first call the ``domain``
852 argument is present so the actual anomaly data can be created.
853 In subsequent calls this may not be true.
854 :note: method needs to be overwritten
855 """
856 raise NotImplementedError()
857
858 class SyntheticFeatureData(SyntheticDataBase):
859 """
860 Uses a list of `SourceFeature` objects to define synthetic anomaly data.
861 """
862 def __init__(self, data_type,
863 features,
864 DIM=2,
865 number_of_elements=10,
866 length=1*U.km,
867 B_b=None,
868 data_offset=0,
869 full_knowledge=False):
870 """
871 :param data_type: data type indicator
872 :type data_type: `DataSource.GRAVITY`, `DataSource.MAGNETIC`
873 :param features: list of features. It is recommended that the features
874 are located entirely below the surface.
875 :type features: ``list`` of `SourceFeature`
876 :param DIM: spatial dimensionality
877 :type DIM: ``int`` (2 or 3)
878 :param number_of_elements: lateral number of elements in the region
879 where data are collected
880 :type number_of_elements: ``int``
881 :param length: lateral extent of the region where data are collected
882 :type length: ``float``
883 :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude]. Only used for magnetic data.
884 :type B_b: ``list`` of ``Scalar``
885 :param data_offset: offset of the data collection region from the surface
886 :type data_offset: ``float``
887 :param full_knowledge: if ``True`` data are collected from the entire subsurface region. This is mainly for testing.
888 :type full_knowledge: ``Bool``
889 """
890 super(SyntheticFeatureData,self).__init__(
891 data_type=data_type, DIM=DIM,
892 number_of_elements=number_of_elements,
893 length=length, B_b=B_b,
894 data_offset=data_offset,
895 full_knowledge=full_knowledge)
896 self._features = features
897
898
899 def getReferenceProperty(self, domain=None):
900 """
901 Returns the reference `Data` object that was used to generate
902 the gravity/susceptibility anomaly data.
903 """
904 if self._reference_data == None:
905 DIM=domain.getDim()
906 x=domain.getX()
907 k=0.
908 for f in self._features:
909 m=f.getMask(x)
910 k = k * (1-m) + f.getValue(x) * m
911 self._reference_data= k
912 return self._reference_data
913
914 class SyntheticData(SyntheticDataBase):
915 """
916 Defines synthetic gravity/magnetic data based on harmonic property anomaly
917
918 rho = amplitude * sin(n_depth * pi /depth * (z+depth_offset)) * sin(n_length * pi /length * (x - shift) )
919
920 for all x and z<=0. For z>0 rho = 0.
921 """
922 def __init__(self, data_type,
923 n_length=1,
924 n_depth=1,
925 depth_offset=0.,
926 depth=None,
927 amplitude=None,
928 DIM=2,
929 number_of_elements=10,
930 length=1*U.km,
931 B_b=None,
932 data_offset=0,
933 full_knowledge=False,
934 s=0.):
935 """
936 :param data_type: data type indicator
937 :type data_type: `DataSource.GRAVITY`, `DataSource.MAGNETIC`
938 :param n_length: number of oscillations in the anomaly data within the
939 observation region
940 :type n_length: ``int``
941 :param n_depth: number of oscillations in the anomaly data below surface
942 :type n_depth: ``int``
943 :param depth_offset: vertical offset of the data
944 :type depth_offset: ``float``
945 :param depth: vertical extent in the anomaly data. If not present the
946 depth of the domain is used.
947 :type depth: ``float``
948 :param amplitude: data amplitude. Default value is 200 U.kg/U.m**3 for
949 gravity and 0.1 for magnetic data.
950 :param DIM: spatial dimensionality
951 :type DIM: ``int`` (2 or 3)
952 :param number_of_elements: lateral number of elements in the region
953 where data are collected
954 :type number_of_elements: ``int``
955 :param length: lateral extent of the region where data are collected
956 :type length: ``float``
957 :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude].
958 Only used for magnetic data.
959 :type B_b: ``list`` of ``Scalar``
960 :param data_offset: offset of the data collection region from the surface
961 :type data_offset: ``float``
962 :param full_knowledge: if ``True`` data are collected from the entire
963 subsurface region. This is mainly for testing.
964 :type full_knowledge: ``Bool``
965 """
966 super(SyntheticData,self).__init__(
967 data_type=data_type, DIM=DIM,
968 number_of_elements=number_of_elements,
969 length=length, B_b=B_b,
970 data_offset=data_offset,
971 full_knowledge=full_knowledge)
972 self.__n_length = n_length
973 self.__n_depth = n_depth
974 self.depth = depth
975 self.depth_offset=depth_offset
976 if amplitude == None:
977 if data_type == DataSource.GRAVITY:
978 amplitude = 200 *U.kg/U.m**3
979 else:
980 amplitude = 0.1
981 self.__amplitude = amplitude
982 self.__s=s
983
984 def getReferenceProperty(self, domain=None):
985 """
986 Returns the reference `Data` object that was used to generate
987 the gravity/susceptibility anomaly data.
988 """
989 if self._reference_data is None:
990 DIM=domain.getDim()
991 x=domain.getX()
992 # set the reference data
993 z=x[DIM-1]
994 dd=self.depth
995 if dd is None: dd=inf(z)
996 z2=(z+self.depth_offset)/(self.depth_offset-dd)
997 k=sin(self.__n_depth * np.pi * z2) * whereNonNegative(z2) * whereNonPositive(z2-1.) * self.__amplitude
998 for i in range(DIM-1):
999 x_i=x[i]
1000 min_x=inf(x_i)
1001 max_x=sup(x_i)
1002 k*= sin(self.__n_length*np.pi*(x_i-min_x-self.__s)/(max_x-min_x))
1003 self._reference_data= k
1004 return self._reference_data
1005
1006 ##############################################################################
1007 class NumpyData(DataSource):
1008 """
1009 """
1010
1011 def __init__(self, data_type, data, error=1., length=1.*U.km, null_value=-1.):
1012 """
1013 A data source that uses survey data from a ``numpy`` object or list
1014 instead of a file.
1015 The dimensionality is inferred from the shape of ``data`` (1- and
1016 2-dimensional data is supported). The data origin is assumed to be
1017 at the coordinate origin.
1018
1019 :param data_type: data type indicator
1020 :type data_type: `DataSource.GRAVITY`, `DataSource.MAGNETIC`
1021 :param data: the survey data array. Note that for a cartesian coordinate
1022 system the shape of the data is considered to be
1023 (nz,ny,nx).
1024 :type data: ``numpy.array`` or ``list``
1025 :param error: constant value to use for the data uncertainties or a
1026 numpy object with uncertainties for every data point.
1027 :type error: ``float`` or ``list`` or ``ndarray``
1028 :param length: side length(s) of the data slice/volume. This can be
1029 a scalar to indicate constant length in all dimensions
1030 or an array/list of values in each coordinate dimension.
1031 :type length: ``float`` or ``list`` or ``ndarray``
1032 :param null_value: value that is used in the undefined regions of the
1033 survey data object.
1034 :type null_value: ``float``
1035 """
1036 super(NumpyData, self).__init__()
1037 if not data_type in [self.GRAVITY, self.MAGNETIC]:
1038 raise ValueError("Invalid value for data_type parameter")
1039 self.__data_type = data_type
1040 self.__data = np.asarray(data, dtype=np.float32)
1041 DIM = len(self.__data.shape)
1042 if DIM not in (1,2):
1043 raise ValueError("NumpyData requires 1- or 2-dimensional data")
1044 self.__error = np.asarray(error, dtype=np.float32)
1045 if len(self.__error.shape) > 0 and \
1046 self.__error.shape != self.__data.shape:
1047 raise ValueError("error argument must be scalar or match the shape of the data")
1048 if isinstance(length,float) or isinstance(length,int):
1049 length = [float(length)] * DIM
1050 else:
1051 length=np.asarray(length, dtype=float)
1052 if len(length.shape) != 1 or length.shape[0] != DIM:
1053 raise ValueError("length must be scalar or an array with %d values"%DIM)
1054 length=length.tolist()
1055
1056 self.__null_value = null_value
1057 self.__nPts = list(reversed(self.__data.shape))
1058 self.__delta = [length[i]/self.__nPts[i] for i in range(DIM)]
1059 self.__origin = [0.] * DIM # this could be an optional argument
1060
1061 def getDataExtents(self):
1062 return (self.__origin, self.__nPts, self.__delta)
1063
1064 def getDataType(self):
1065 return self.__data_type
1066
1067 def getSurveyData(self, domain, origin, NE, spacing):
1068 FS = ReducedFunction(domain)
1069 nValues = self.__nPts
1070 dataDIM = len(nValues)
1071 # determine base location of this dataset within the domain
1072 first=[int((self.__origin[i]-origin[i])/spacing[i]) for i in range(dataDIM)]
1073 # determine the resolution difference between domain and data.
1074 # If domain has twice the resolution we can double up the data etc.
1075 multiplier=[int(round(self.__delta[i]/spacing[i])) for i in range(dataDIM)]
1076 if domain.getDim() > dataDIM:
1077 first.append(int((-origin[-1])/spacing[-1]))
1078 multiplier=multiplier+[1]
1079 nValues=nValues+[1]
1080
1081 _handle, numpyfile = tempfile.mkstemp()
1082 os.close(_handle)
1083 self.__data.tofile(numpyfile)
1084
1085 data = ripleycpp._readBinaryGrid(numpyfile, FS, first, nValues,
1086 multiplier, (), self.__null_value)
1087 if len(self.__error.shape) > 0:
1088 self.__error.tofile(numpyfile)
1089 sigma = ripleycpp._readBinaryGrid(numpyfile, FS, first, nValues,
1090 multiplier, (), 0.)
1091 else:
1092 sigma = self.__error.item() * whereNonZero(data-self.__null_value)
1093
1094 os.unlink(numpyfile)
1095
1096 return data, sigma
1097
1098 def getUtmZone(self):
1099 """
1100 returns a dummy UTM zone since this class does not use real coordinate
1101 values.
1102 """
1103 return 0
1104

  ViewVC Help
Powered by ViewVC 1.1.26