Package eoxserver :: Package services :: Package ows :: Package wcst :: Module gdalGeoTiff
Source Code for Module eoxserver.services.ows.wcst.gdalGeoTiff

  1  #!/usr/bin/env python  
  2  #----------------------------------------------------------------------- 
  3  # $Id: gdalGeoTiff.py 2288 2013-03-08 15:56:55Z meissls $ 
  4  # 
  5  # Description:  
  6  # 
  7  # GDAL based GeoTIFF utilities 
  8  #  
  9  #------------------------------------------------------------------------------- 
 10  # 
 11  # Project: EOxServer <http://eoxserver.org> 
 12  # Authors: Martin Paces <martin.paces@iguassu.cz> 
 13  # 
 14  #------------------------------------------------------------------------------- 
 15  # Copyright (C) 2011 Iguassu Software Systems, a.s  
 16  # 
 17  # Permission is hereby granted, free of charge, to any person obtaining a copy 
 18  # of this software and associated documentation files (the "Software"), to deal 
 19  # in the Software without restriction, including without limitation the rights 
 20  # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell  
 21  # copies of the Software, and to permit persons to whom the Software is  
 22  # furnished to do so, subject to the following conditions: 
 23  # 
 24  # The above copyright notice and this permission notice shall be included in all 
 25  # copies of this Software or works derived from this Software. 
 26  # 
 27  # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 
 28  # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 
 29  # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 
 30  # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 
 31  # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 
 32  # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 
 33  # THE SOFTWARE. 
 34  #------------------------------------------------------------------------------- 
 35   
 36  import sys 
 37  import os.path  
 38  from numpy import linalg  
 39  from numpy import array 
 40  from numpy import dot  
 41  from numpy import argmin, argmax  
 42   
 43  from qhull import chull2D_qhull 
 44   
 45  from eoxserver.contrib import osr, gdal 
 46   
 47  #------------------------------------------------------------------------------- 
 48   
 49 -class GeoTransform(object) :  
 50   
 51 -    def __init__( self , gt ) :  
 52           
 53          self.xy = array( (gt[0],gt[3]) , 'float64' )  
 54          self.af = array( ((gt[2],gt[1]),(gt[5],gt[4])) , 'float64' )  
 55          self.ab = linalg.inv( self.af )  
 56   
 57 -    def rc2xy( self , rc ) :  
 58          rc = array( rc , 'float64' )  
 59          xy = self.xy + dot( self.af , rc )   
 60          return xy  
 61   
 62 -    def xy2rc( self , xy ) :  
 63          xy = array( xy , 'float64' )  
 64          rc = dot( self.ab , xy - self.xy )   
 65          return rc  
 66   
 67  #------------------------------------------------------------------------------- 
 68   
 69 -class OSRTransform( object ) :  
 70   
 71 -    def __init__( self , src , dst ) :  
 72   
 73          self.fwd = osr.CoordinateTransformation( src , dst ) 
 74          self.bwd = osr.CoordinateTransformation( dst , src ) 
 75   
 76 -    def src2dst( self , xy ) : 
 77          return array( self.fwd.TransformPoint( xy[0] , xy[1] ) , 'float64' )  
 78   
 79 -    def dst2src( self , xy ) : 
 80          return array( self.bwd.TransformPoint( xy[0] , xy[1] ) , 'float64' )  
 81   
 82  #------------------------------------------------------------------------------- 
 83   
 84 -class GDalInfo( object ) :  
 85      """ this class holds the GDAL information of the rester file """ 
 86   
 87 -    def __init__( self , fname ) :  
 88   
 89          # get the dataset and its params 
 90          gdal.AllRegister() 
 91          ds = gdal.Open(fname) 
 92       
 93          self.fileName = fname  
 94          self.driverName =  ds.GetDriver().LongName 
 95          self.size       = (ds.RasterYSize,ds.RasterXSize,ds.RasterCount)  
 96          self.GCPCount   = ds.GetGCPCount()  
 97          self.GCPProjection = ds.GetGCPProjection() 
 98          self.Projection = ds.GetProjection() 
 99          self.ProjectionRef = ds.GetProjectionRef()  
100          self.GeoTransform = ds.GetGeoTransform()  
101          self.GCP         = ds.GetGCPs()  
102          self.isRectified  = bool( self.Projection )  
103          self.isReferenceable = bool( self.GCPProjection) and ( self.GCPCount > 0 )  
104   
105   
106 -    def __str__( self ) :  
107   
108          a = []  
109          a.append( "GDalInfo:" )  
110          a.append( "\tFile:         \t %s" % str(self.fileName) )  
111          a.append( "\tDriver:       \t %s" % str(self.driverName) )  
112          a.append( "\tSize:         \t %s" % str( self.size ) ) 
113          a.append( "\tNo.GCP        \t %s" % str( self.GCPCount ) ) 
114          a.append( "\tProjectionRef:\t %s" % str( self.ProjectionRef  ) ) 
115          a.append( "\tProjection:   \t %s" % str( self.Projection ) ) 
116          a.append( "\tGeoTransform: \t %s" % str( self.GeoTransform ) ) 
117          a.append( "\tGCP Projection:\t %s" % str( self.GCPProjection ) ) 
118          a.append( "" )  
119   
120          return "\n".join( a )  
121   
122  #------------------------------------------------------------------------------- 
123   
124 -def getFootprint( info , **kvarg ) :  
125      """ extract geotiff image footprint  
126       
127          info            - [instance of GDalInfo class] image descriptor   
128          numberOfPoints  - [integer] number of points per single image size  
129                            (20 by default, applicable to Rectified images only)  
130          delimiter       - [string] separator delimiting the footprint coordinates (<space> by default)  
131          repeatFirst     - [boolean] if True the first point is repeated as last one to close the loop 
132      """ 
133   
134      if info.isRectified :  
135          return getFootprintRect( info , **kvarg ) 
136      elif info.isReferenceable :  
137          return getFootprintRef( info , **kvarg ) 
138   
139      return None  
140   
141   
142 -def getFootprintRect( info , numberOfPoints = 20, delimiter = " " , repeatFirst = False ) :  
143      """ extract geotiff image footprint of a Rectified image   
144       
145          info            - [instance of GDalInfo class] image descriptor   
146          numberOfPoints  - [integer] number of points per single image size  
147                            (20 by default)  
148          delimiter       - [string] separator delimiting the footprint coordinates (<space> by default)  
149          repeatFirst     - [boolean] if True the first point is repeated as last one to close the loop 
150      """ 
151   
152      # ----------------------------- 
153      # projections conversion 
154   
155      size = ( info.size[0] , info.size[1] )  
156   
157      gt = GeoTransform( info.GeoTransform )  
158   
159      cr_src = osr.SpatialReference()  
160      cr_dst = osr.SpatialReference()  
161   
162      cr_src.ImportFromWkt( info.Projection )  
163      cr_dst.SetWellKnownGeogCS( "WGS84" )  
164   
165      trn = OSRTransform(cr_src,cr_dst)  
166   
167      # ----------------------------- 
168      # extract footprint (clockwise) 
169   
170      foot = []  
171   
172      # top (L2R) 
173      for i  in xrange( 0 , numberOfPoints ) : foot.append( tuple( trn.src2dst(gt.rc2xy((0,size[1]*float(i)/float(numberOfPoints))))[:2] ) )  
174      # right (T2B)  
175      for i  in xrange( 0 , numberOfPoints ) : foot.append( tuple( trn.src2dst(gt.rc2xy((size[0]*float(i)/float(numberOfPoints),size[1])))[:2] ) )  
176      # bottom (R2L)  
177      for i  in xrange( 0 , numberOfPoints ) : foot.append( tuple( trn.src2dst(gt.rc2xy((size[0],size[1]*float(numberOfPoints-i)/float(numberOfPoints))))[:2] ) )  
178      # left (B2T)  
179      for i  in xrange( 0 , numberOfPoints ) : foot.append( tuple( trn.src2dst(gt.rc2xy((size[0]*float(numberOfPoints-i)/float(numberOfPoints),0)))[:2] ) )  
180   
181      if ( repeatFirst ) : foot.append( foot[0] )  
182       
183      # things to be done  
184      # 1) POLAR PROJECTIONS  
185      # 2) NON-POLAR FOOTPRINT CROSSING +/-180 MERIDIAN  
186   
187      return delimiter.join(map( lambda x : "%.6f%s%.6f"%(x[0],delimiter,x[1]) , foot )) 
188   
189  #------------------------------------------------------------------------------- 
190   
191 -def getFootprintRef( info , numberOfPoints = 20, delimiter = " " , repeatFirst = False ) :  
192      """ extract geotiff image footprint of a Referenceableimage   
193       
194          info            - [instance of GDalInfo class] image descriptor   
195          numberOfPoints  - [integer] dummy placeholder - not used (20 by default)  
196          delimiter       - [string] separator delimiting the footprint coordinates (<space> by default)  
197          repeatFirst     - [boolean] if True the first point is repeated as last one to close the loop 
198      """ 
199      # NOTE: numberOfPoints parameter is ignored  
200   
201      D = array( map( lambda gcp : ( gcp.GCPLine , gcp.GCPPixel ) , info.GCP ) )  
202   
203      H = chull2D_qhull( D , eps = 1e-3 , clock_wise = True )  
204   
205      if ( repeatFirst ) : H.append( H[0] )  
206   
207      # things to be done  
208      # 1) POLAR PROJECTIONS  
209      # 2) NON-POLAR FOOTPRINT CROSSING +/-180 MERIDIAN  
210   
211      return delimiter.join(map( lambda i : "%.6f%s%.6f"%(info.GCP[i].GCPX,delimiter,info.GCP[i].GCPY) , H )) 
212   
213  #------------------------------------------------------------------------------- 
214   
215  # simple test  
216  if __name__ == "__main__" :  
217           
218      fname = os.path.abspath( sys.argv[1] )  
219   
220      info = GDalInfo( fname )  
221   
222      #print info  
223      #print "WGS84 Footprint: "  
224      print getFootprint( info ) 
225