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

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

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

revision 3947 by caltinay, Wed Aug 22 23:19:10 2012 UTC revision 3956 by caltinay, Wed Sep 5 04:57:43 2012 UTC
# Line 23  __all__ = ['DataSource','UBCDataSource', Line 23  __all__ = ['DataSource','UBCDataSource',
23    
24  import logging  import logging
25  import numpy as np  import numpy as np
26    import pyproj
27  from esys.escript import *  from esys.escript import *
28  from esys.escript.linearPDEs import *  from esys.escript.linearPDEs import *
29  import esys.escript.unitsSI as U  import esys.escript.unitsSI as U
# Line 32  try: Line 33  try:
33  except:  except:
34      pass      pass
35    
36    def LatLonToUTM(lon, lat, wkt_string):
37        zone=int(np.median((np.floor((np.array(lon) + 180)/6) + 1) % 60))
38        try:
39            import osgeo.osr
40            srs = osgeo.osr.SpatialReference()
41            srs.ImportFromWkt(wkt_string)
42            p_src = pyproj.Proj(srs.ExportToProj4())
43        except:
44            p_src = pyproj.Proj('+proj=longlat +ellps=clrk66 +no_defs')
45        p_dest = pyproj.Proj(proj='utm', zone=zone) # ellps?
46        x,y=pyproj.transform(p_src, p_dest, lon, lat)
47        return x,y
48    
49  class DataSource(object):  class DataSource(object):
50      """      """
# Line 234  class UBCDataSource(DataSource): Line 247  class UBCDataSource(DataSource):
247    
248  ##############################################################################  ##############################################################################
249  class NetCDFDataSource(DataSource):  class NetCDFDataSource(DataSource):
250      def __init__(self, domainclass, gravfile, topofile=None, vertical_extents=(-40000,10000,26), alt_of_data=0.):      def __init__(self, domainclass, gravfile, topofile=None, vertical_extents=(-40000,10000,26), alt_of_data=1.):
251          """          """
252          vertical_extents - (alt_min, alt_max, num_elements)          vertical_extents - (alt_min, alt_max, num_elements)
253          alt_of_data - altitude of measurements          alt_of_data - altitude of measurements
# Line 278  class NetCDFDataSource(DataSource): Line 291  class NetCDFDataSource(DataSource):
291              raise RuntimeError("Could not determine extents of data")              raise RuntimeError("Could not determine extents of data")
292    
293          # find longitude and latitude variables          # find longitude and latitude variables
         variables=f.variables  
294          lon_name=None          lon_name=None
295          for n in ['lon','longitude']:          for n in ['lon','longitude']:
296              if n in variables:              if n in f.variables:
297                  lon_name=n                  lon_name=n
298                  variables.pop(n)                  longitude=f.variables.pop(n)
299                  break                  break
300          if lon_name is None:          if lon_name is None:
301              raise RuntimeError("Could not determine longitude variable")              raise RuntimeError("Could not determine longitude variable")
302          lat_name=None          lat_name=None
303          for n in ['lat','latitude']:          for n in ['lat','latitude']:
304              if n in variables:              if n in f.variables:
305                  lat_name=n                  lat_name=n
306                  variables.pop(n)                  latitude=f.variables.pop(n)
307                  break                  break
308          if lat_name is None:          if lat_name is None:
309              raise RuntimeError("Could not determine latitude variable")              raise RuntimeError("Could not determine latitude variable")
310    
311          # try to figure out gravity variable name          # try to figure out gravity variable name
312          grav_name=None          grav_name=None
313          if len(variables)==1:          if len(f.variables)==1:
314              grav_name=variables.keys()[0]              grav_name=f.variables.keys()[0]
315          else:          else:
316              for n in variables.keys():              for n in f.variables.keys():
317                  dims=variables[n].dimensions                  dims=f.variables[n].dimensions
318                  if (lat_name in dims) and (lon_name in dims):                  if (lat_name in dims) and (lon_name in dims):
319                      grav_name=n                      grav_name=n
320                      break                      break
321          if grav_name is None:          if grav_name is None:
322              raise RuntimeError("Could not determine gravity variable")              raise RuntimeError("Could not determine gravity variable")
323    
324            # see if there is a wkt string to convert coordinates
325            try:
326                wkt_string=f.variables[grav_name].esri_pe_string
327            except:
328                wkt_string=None
329    
330          origin=[0.,0.,ve[0]]          origin=[0.,0.,ve[0]]
331          lengths=[1.,1.,ve[1]-ve[0]]          lengths=[100000.,100000.,ve[1]-ve[0]]
332    
333          try:          try:
334              lon_range=f.variables[lon_name].actual_range              lon_range=longitude.actual_range
335              lat_range=f.variables[lat_name].actual_range              lat_range=latitude.actual_range
336              #origin[:2]=[lon_range[0],lat_range[0]]              lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)
337                origin[:2]=lon_range[0],lat_range[0]
338              lengths[:2]=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]              lengths[:2]=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]
339          except:          except:
340              try:              try:
341                  #origin[:2]=[f.geospatial_lon_min, f.geospatial_lat_min]                  lon_range=[f.geospatial_lon_min,f.geospatial_lon_max]
342                  lengths[:2]=[f.geospatial_lon_max-origin[0], f.geospatial_lat_max-origin[1]]                  lat_range=[f.geospatial_lat_min,f.geospatial_lat_max]
343                    lon_range,lat_range=LatLonToUTM(lon_range, lat_range, wkt_string)
344                    origin[:2]=lon_range[0],lat_range[0]
345                    lengths[:2]=[lon_range[1]-lon_range[0], lat_range[1]-lat_range[0]]
346              except:              except:
347                  pass                  pass
348    
         # hack  
         # grid spacing ~800m  
         SPACING_X=SPACING_Y=800.  
         lengths[0]=(NX-1)*SPACING_X  
         lengths[1]=(NY-1)*SPACING_Y  
349          f.close()          f.close()
350    
351          self._numDataPoints=[NX, NY, ve[2]]          self._numDataPoints=[NX, NY, ve[2]]
352          self._origin=origin          self._origin=origin
353          self._meshlen=lengths          self._spacing=[np.round(lengths[i]/(self._numDataPoints[i]-1)) for i in xrange(3)]
354          self._spacing=[lengths[i]/(self._numDataPoints[i]-1) for i in xrange(3)]          self._meshlen=[self._numDataPoints[i]*self._spacing[i] for i in xrange(3)]
355            self._wkt_string=wkt_string
356          self._lon=lon_name          self._lon=lon_name
357          self._lat=lat_name          self._lat=lat_name
358          self._grv=grav_name          self._grv=grav_name
# Line 352  class NetCDFDataSource(DataSource): Line 372  class NetCDFDataSource(DataSource):
372              dom=self._domainclass(*self.NE, l0=lo[0], l1=lo[1], l2=lo[2])              dom=self._domainclass(*self.NE, l0=lo[0], l1=lo[1], l2=lo[2])
373              # ripley may adjust NE and length              # ripley may adjust NE and length
374              self._meshlen=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in xrange(3)]              self._meshlen=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in xrange(3)]
375              self.NE=[self._meshlen[i]/(self._spacing[i]) for i in xrange(3)]              self.NE=[self._meshlen[i]/self._spacing[i] for i in xrange(3)]
376              x=dom.getX()-[self._origin[i]+NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]              x=dom.getX()-[self._origin[i]+NEdiff[i]/2.*self._spacing[i] for i in xrange(3)]
377              mask=wherePositive(dom.getX()[2])              mask=wherePositive(dom.getX()[2])
378    
# Line 404  class NetCDFDataSource(DataSource): Line 424  class NetCDFDataSource(DataSource):
424    
425      def _readGravity(self):      def _readGravity(self):
426          f=netcdf_file(self._gravfile, 'r')          f=netcdf_file(self._gravfile, 'r')
427          #lon=f.variables[self._lon][:]          lon=f.variables[self._lon][:]
428          #lat=f.variables[self._lat][:]          lat=f.variables[self._lat][:]
429            lon,lat=np.meshgrid(lon,lat)
430            lon,lat=LatLonToUTM(lon, lat, self._wkt_string)
431          grav=f.variables[self._grv][:]          grav=f.variables[self._grv][:]
432          f.close()          f.close()
         lon=np.linspace(0, (grav.shape[1]-1)*self._spacing[0], num=grav.shape[1])  
         lat=np.linspace(0, (grav.shape[0]-1)*self._spacing[1], num=grav.shape[0])  
         lon,lat=np.meshgrid(lon,lat)  
433          lon=lon.flatten()          lon=lon.flatten()
434          lat=lat.flatten()          lat=lat.flatten()
435          grav=grav.flatten()          grav=grav.flatten()
436          alt=self._altOfData*np.ones(grav.shape)          alt=self._altOfData*np.ones(grav.shape)
437            # error value is an assumption
438          try:          try:
439              missing=grav.missing_value              missing=grav.missing_value
440              err=np.where(grav==missing, 1.0, 0.0)              err=np.where(grav==missing, 20.0, 0.0)
441          except:          except:
442              err=np.ones(lon.shape)              err=20.0*np.ones(lon.shape)
443            # convert units
444          err=1e-6*err          err=1e-6*err
445          grav=1e-6*grav          grav=1e-6*grav
446          gravdata=np.column_stack((lon,lat,alt,grav,err))          gravdata=np.column_stack((lon,lat,alt,grav,err))

Legend:
Removed from v.3947  
changed lines
  Added in v.3956

  ViewVC Help
Powered by ViewVC 1.1.26