/[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 3958 by caltinay, Fri Sep 7 02:56:42 2012 UTC revision 3959 by caltinay, Fri Sep 7 04:48:07 2012 UTC
# Line 59  def LatLonToUTM(lon, lat, wkt_string=Non Line 59  def LatLonToUTM(lon, lat, wkt_string=Non
59          p_src = pyproj.Proj(srs.ExportToProj4())          p_src = pyproj.Proj(srs.ExportToProj4())
60      except:      except:
61          p_src = pyproj.Proj('+proj=longlat +ellps=clrk66 +no_defs')          p_src = pyproj.Proj('+proj=longlat +ellps=clrk66 +no_defs')
62      p_dest = pyproj.Proj(proj='utm', zone=zone) # ellps?      # we assume southern hemisphere here
63        p_dest = pyproj.Proj('+proj=utm +zone=%d +south +units=m'%zone)
64      x,y=pyproj.transform(p_src, p_dest, lon, lat)      x,y=pyproj.transform(p_src, p_dest, lon, lat)
65      return x,y      return x,y
66    
# Line 114  class DataSource(object): Line 115  class DataSource(object):
115          for entry in data:          for entry in data:
116              index=()              index=()
117              for i in xrange(dim):              for i in xrange(dim):
118                  index=((entry[i]-origin[i])/spacing[i],)+index                  index=(int((entry[i]-origin[i])/spacing[i]),)+index
119              for i in xrange(arrays.shape[0]):              for i in xrange(arrays.shape[0]):
120                  arrays[i][index]=entry[dim+i]                  arrays[i][index]=entry[dim+i]
121          dom=self.getDomain()          dom=self.getDomain()
# Line 547  class ERSDataSource(DataSource): Line 548  class ERSDataSource(DataSource):
548              raise RuntimeError("Could not determine cell dimensions")              raise RuntimeError("Could not determine cell dimensions")
549    
550          try:          try:
551              originX = float(md_dict['RasterInfo.RegistrationCoord.Eastings'])              if md_dict['CoordinateSpace.CoordinateType']=='EN':
552              originY = float(md_dict['RasterInfo.RegistrationCoord.Northings'])                  originX = float(md_dict['RasterInfo.RegistrationCoord.Eastings'])
553                    originY = float(md_dict['RasterInfo.RegistrationCoord.Northings'])
554                elif md_dict['CoordinateSpace.CoordinateType']=='LL':
555                    originX = float(md_dict['RasterInfo.RegistrationCoord.Longitude'])
556                    originY = float(md_dict['RasterInfo.RegistrationCoord.Latitude'])
557                else:
558                    raise RuntimeError("Unknown CoordinateType")
559          except:          except:
560              self.logger.warn("Could not determine coordinate origin")              self.logger.warn("Could not determine coordinate origin. Setting to (0.0, 0.0)")
561              originX,originY = 0.0, 0.0              originX,originY = 0.0, 0.0
562    
563          # test data sets have origin in top-left corner so y runs top-down          if 'GEODETIC' in md_dict['CoordinateSpace.Projection']:
564                # it appears we have lat/lon coordinates so need to convert
565                # origin and spacing. Try using gdal to get the wkt if available:
566                try:
567                    from osgeo import gdal
568                    ds=gdal.Open(self._headerfile)
569                    wkt=ds.GetProjection()
570                except:
571                    wkt='GEOGCS["GEOCENTRIC DATUM of AUSTRALIA",DATUM["GDA94",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'
572                    self.logger.warn('GDAL not available, assuming GDA94 data')
573                originX_UTM,originY_UTM=LatLonToUTM(originX, originY, wkt)
574                op1X,op1Y=LatLonToUTM(originX+spacingX, originY+spacingY, wkt)
575                # we are rounding to avoid interpolateOnDomain issues
576                spacingX=np.round(op1X-originX_UTM)
577                spacingY=np.round(op1Y-originY_UTM)
578                originX=np.round(1000*originX_UTM)/1000.
579                originY=np.round(1000*originY_UTM)/1000.
580    
581            # data sets have origin in top-left corner so y runs top-down
582          self._dataorigin=[originX, originY]          self._dataorigin=[originX, originY]
583          originY-=(NY-1)*spacingY          originY-=(NY-1)*spacingY
584          self._maskval=maskval          self._maskval=maskval
# Line 623  class ERSDataSource(DataSource): Line 648  class ERSDataSource(DataSource):
648    
649          NE=[self._numDataPoints[i] for i in xrange(2)]          NE=[self._numDataPoints[i] for i in xrange(2)]
650          x=self._dataorigin[0]+np.arange(start=0, stop=NE[0]*self._spacing[0], step=self._spacing[0])          x=self._dataorigin[0]+np.arange(start=0, stop=NE[0]*self._spacing[0], step=self._spacing[0])
651          # test data sets have origin in top-left corner so y runs top-down          # data sets have origin in top-left corner so y runs top-down
652          y=self._dataorigin[1]-np.arange(start=0, stop=NE[1]*self._spacing[1], step=self._spacing[1])          y=self._dataorigin[1]-np.arange(start=0, stop=NE[1]*self._spacing[1], step=self._spacing[1])
653          x,y=np.meshgrid(x,y)          x,y=np.meshgrid(x,y)
654          x,y=x.flatten(),y.flatten()          x,y=x.flatten(),y.flatten()

Legend:
Removed from v.3958  
changed lines
  Added in v.3959

  ViewVC Help
Powered by ViewVC 1.1.26