/[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 4561 - (show annotations)
Wed Dec 4 05:19:32 2013 UTC (5 years, 10 months ago) by caltinay
File MIME type: text/x-python
File size: 49412 byte(s)
Removed spurious print and some minor changes in inversion example.

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

  ViewVC Help
Powered by ViewVC 1.1.26