1
2
3
4
5
6
7
8 from sys import version_info
9 if version_info >= (2,6,0):
11 from os.path import dirname
12 import imp
13 fp = None
14 try:
15 fp, pathname, description = imp.find_module('_gdal_array', [dirname(__file__)])
16 except ImportError:
17 import _gdal_array
18 return _gdal_array
19 if fp is not None:
20 try:
21 _mod = imp.load_module('_gdal_array', fp, pathname, description)
22 finally:
23 fp.close()
24 return _mod
25 _gdal_array = swig_import_helper()
26 del swig_import_helper
27 else:
28 import _gdal_array
29 del version_info
30 try:
31 _swig_property = property
32 except NameError:
33 pass
35 if (name == "thisown"): return self.this.own(value)
36 if (name == "this"):
37 if type(value).__name__ == 'SwigPyObject':
38 self.__dict__[name] = value
39 return
40 method = class_type.__swig_setmethods__.get(name,None)
41 if method: return method(self,value)
42 if (not static) or hasattr(self,name):
43 self.__dict__[name] = value
44 else:
45 raise AttributeError("You cannot add attributes to %s" % self)
46
49
51 if (name == "thisown"): return self.this.own()
52 method = class_type.__swig_getmethods__.get(name,None)
53 if method: return method(self)
54 raise AttributeError(name)
55
57 try: strthis = "proxy of " + self.this.__repr__()
58 except: strthis = ""
59 return "<%s.%s; %s >" % (self.__class__.__module__, self.__class__.__name__, strthis,)
60
61 try:
62 _object = object
63 _newclass = 1
64 except AttributeError:
66 _newclass = 0
67
68
69
71 """GetArrayFilename(PyArrayObject psArray) -> retStringAndCPLFree"""
72 return _gdal_array.GetArrayFilename(*args)
73
75 """
76 BandRasterIONumPy(Band band, int bWrite, int xoff, int yoff, int xsize,
77 int ysize, PyArrayObject psArray, int buf_type) -> CPLErr
78 """
79 return _gdal_array.BandRasterIONumPy(*args, **kwargs)
80 import numpy
81 import _gdal_array
82
83 import gdalconst
84 import gdal
85 gdal.AllRegister()
86
87 codes = { gdalconst.GDT_Byte : numpy.uint8,
88 gdalconst.GDT_UInt16 : numpy.uint16,
89 gdalconst.GDT_Int16 : numpy.int16,
90 gdalconst.GDT_UInt32 : numpy.uint32,
91 gdalconst.GDT_Int32 : numpy.int32,
92 gdalconst.GDT_Float32 : numpy.float32,
93 gdalconst.GDT_Float64 : numpy.float64,
94 gdalconst.GDT_CInt16 : numpy.complex64,
95 gdalconst.GDT_CInt32 : numpy.complex64,
96 gdalconst.GDT_CFloat32 : numpy.complex64,
97 gdalconst.GDT_CFloat64 : numpy.complex128
98 }
99
101
102 ds = gdal.Open( GetArrayFilename(array) )
103
104 if ds is not None and prototype_ds is not None:
105 if type(prototype_ds).__name__ == 'str':
106 prototype_ds = gdal.Open( prototype_ds )
107 if prototype_ds is not None:
108 CopyDatasetInfo( prototype_ds, ds )
109
110 return ds
111
112
114 if isinstance(code, type):
115
116
117 if code == numpy.int8:
118 return gdalconst.GDT_Byte
119 if code == numpy.complex64:
120 return gdalconst.GDT_CFloat32
121
122 for key, value in codes.items():
123 if value == code:
124 return key
125 return None
126 else:
127 try:
128 return codes[code]
129 except KeyError:
130 return None
131
133 if not isinstance(numeric_type, type):
134 raise TypeError("Input must be a type")
135 return flip_code(numeric_type)
136
139
140 -def LoadFile( filename, xoff=0, yoff=0, xsize=None, ysize=None ):
146
147 -def SaveArray( src_array, filename, format = "GTiff", prototype = None ):
148 driver = gdal.GetDriverByName( format )
149 if driver is None:
150 raise ValueError("Can't find driver "+format)
151
152 return driver.CreateCopy( filename, OpenArray(src_array,prototype) )
153
155
156 if xsize is None:
157 xsize = ds.RasterXSize
158 if ysize is None:
159 ysize = ds.RasterYSize
160
161 if ds.RasterCount == 1:
162 return BandReadAsArray( ds.GetRasterBand(1), xoff, yoff, xsize, ysize, buf_obj = buf_obj)
163
164 datatype = ds.GetRasterBand(1).DataType
165 for band_index in range(2,ds.RasterCount+1):
166 if datatype != ds.GetRasterBand(band_index).DataType:
167 datatype = gdalconst.GDT_Float32
168
169 typecode = GDALTypeCodeToNumericTypeCode( datatype )
170 if typecode == None:
171 datatype = gdalconst.GDT_Float32
172 typecode = numpy.float32
173
174 if buf_obj is not None:
175 for band_index in range(1,ds.RasterCount+1):
176 BandReadAsArray( ds.GetRasterBand(band_index),
177 xoff, yoff, xsize, ysize, buf_obj = buf_obj[band_index-1])
178 return buf_obj
179
180 array_list = []
181 for band_index in range(1,ds.RasterCount+1):
182 band_array = BandReadAsArray( ds.GetRasterBand(band_index),
183 xoff, yoff, xsize, ysize)
184 array_list.append( numpy.reshape( band_array, [1,ysize,xsize] ) )
185
186 return numpy.concatenate( array_list )
187
188 -def BandReadAsArray( band, xoff = 0, yoff = 0, win_xsize = None, win_ysize = None,
189 buf_xsize=None, buf_ysize=None, buf_obj=None ):
190 """Pure python implementation of reading a chunk of a GDAL file
191 into a numpy array. Used by the gdal.Band.ReadAsArray method."""
192
193 if win_xsize is None:
194 win_xsize = band.XSize
195 if win_ysize is None:
196 win_ysize = band.YSize
197
198 if buf_obj is None:
199 if buf_xsize is None:
200 buf_xsize = win_xsize
201 if buf_ysize is None:
202 buf_ysize = win_ysize
203 else:
204 if len(buf_obj.shape) == 2:
205 shape_buf_xsize = buf_obj.shape[1]
206 shape_buf_ysize = buf_obj.shape[0]
207 else:
208 shape_buf_xsize = buf_obj.shape[2]
209 shape_buf_ysize = buf_obj.shape[1]
210 if buf_xsize is not None and buf_xsize != shape_buf_xsize:
211 raise ValueError('Specified buf_xsize not consistant with array shape')
212 if buf_ysize is not None and buf_ysize != shape_buf_ysize:
213 raise ValueError('Specified buf_ysize not consistant with array shape')
214 buf_xsize = shape_buf_xsize
215 buf_ysize = shape_buf_ysize
216
217 if buf_obj is None:
218 datatype = band.DataType
219 typecode = GDALTypeCodeToNumericTypeCode( datatype )
220 if typecode == None:
221 datatype = gdalconst.GDT_Float32
222 typecode = numpy.float32
223 else:
224 datatype = NumericTypeCodeToGDALTypeCode( typecode )
225
226 if datatype == gdalconst.GDT_Byte and band.GetMetadataItem('PIXELTYPE', 'IMAGE_STRUCTURE') == 'SIGNEDBYTE':
227 typecode = numpy.int8
228 ar = numpy.empty([buf_ysize,buf_xsize], dtype = typecode)
229 if BandRasterIONumPy( band, 0, xoff, yoff, win_xsize, win_ysize,
230 ar, datatype ) != 0:
231 return None
232
233 return ar
234 else:
235 datatype = NumericTypeCodeToGDALTypeCode( buf_obj.dtype.type )
236 if not datatype:
237 raise ValueError("array does not have corresponding GDAL data type")
238
239 if BandRasterIONumPy( band, 0, xoff, yoff, win_xsize, win_ysize,
240 buf_obj, datatype ) != 0:
241 return None
242
243 return buf_obj
244
246 """Pure python implementation of writing a chunk of a GDAL file
247 from a numpy array. Used by the gdal.Band.WriteArray method."""
248
249 if array is None or len(array.shape) != 2:
250 raise ValueError("expected array of dim 2")
251
252 xsize = array.shape[1]
253 ysize = array.shape[0]
254
255 if xsize + xoff > band.XSize or ysize + yoff > band.YSize:
256 raise ValueError("array larger than output file, or offset off edge")
257
258 datatype = NumericTypeCodeToGDALTypeCode( array.dtype.type )
259
260
261
262 if not datatype:
263 gdal.Debug( 'gdal_array', 'force array to float64' )
264 array = array.astype( numpy.float64 )
265 datatype = NumericTypeCodeToGDALTypeCode( array.dtype.type )
266
267 if not datatype:
268 raise ValueError("array does not have corresponding GDAL data type")
269
270 return BandRasterIONumPy( band, 1, xoff, yoff, xsize, ysize,
271 array, datatype )
272
273
275 """
276 Copy georeferencing information and metadata from one dataset to another.
277 src: input dataset
278 dst: output dataset - It can be a ROI -
279 xoff, yoff: dst's offset with respect to src in pixel/line.
280
281 Notes: Destination dataset must have update access. Certain formats
282 do not support creation of geotransforms and/or gcps.
283
284 """
285
286 dst.SetMetadata( src.GetMetadata() )
287
288
289
290
291 gt = src.GetGeoTransform()
292 if gt != (0,1,0,0,0,1):
293 dst.SetProjection( src.GetProjectionRef() )
294
295 if (xoff == 0) and (yoff == 0):
296 dst.SetGeoTransform( gt )
297 else:
298 ngt = [gt[0],gt[1],gt[2],gt[3],gt[4],gt[5]]
299 ngt[0] = gt[0] + xoff*gt[1] + yoff*gt[2];
300 ngt[3] = gt[3] + xoff*gt[4] + yoff*gt[5];
301 dst.SetGeoTransform( ( ngt[0], ngt[1], ngt[2], ngt[3], ngt[4], ngt[5] ) )
302
303
304 elif src.GetGCPCount() > 0:
305
306 if (xoff == 0) and (yoff == 0):
307 dst.SetGCPs( src.GetGCPs(), src.GetGCPProjection() )
308 else:
309 gcps = src.GetGCPs()
310
311 new_gcps = []
312 for gcp in gcps:
313 ngcp = gdal.GCP()
314 ngcp.GCPX = gcp.GCPX
315 ngcp.GCPY = gcp.GCPY
316 ngcp.GCPZ = gcp.GCPZ
317 ngcp.GCPPixel = gcp.GCPPixel - xoff
318 ngcp.GCPLine = gcp.GCPLine - yoff
319 ngcp.Info = gcp.Info
320 ngcp.Id = gcp.Id
321 new_gcps.append(ngcp)
322
323 try:
324 dst.SetGCPs( new_gcps , src.GetGCPProjection() )
325 except:
326 print ("Failed to set GCPs")
327 return
328
329 return
330