我们从Python开源项目中，提取了以下21个代码示例，用于说明如何使用osgeo.gdal.GDT_Float32()。
def write_stats(self): #if not hasattr(self, 'stack_count'): stat_list = ['stack_count', 'stack_mean', 'stack_std', 'stack_min', 'stack_max'] if self.med: stat_list.append('stack_med') stat_list.append('stack_nmad') if any([not hasattr(self, i) for i in stat_list]): self.compute_stats() print("Writing out stats") #Create dummy ds - might want to use vrt here instead driver = gdal.GetDriverByName("MEM") ds = driver.Create('', self.ma_stack.shape[2], self.ma_stack.shape[1], 1, gdal.GDT_Float32) ds.SetGeoTransform(self.gt) ds.SetProjection(self.proj) #Write out with malib, should preserve ma type out_prefix = os.path.splitext(self.stack_fn)[0] iolib.writeGTiff(self.stack_count, out_prefix+'_count.tif', ds) iolib.writeGTiff(self.stack_mean, out_prefix+'_mean.tif', ds) iolib.writeGTiff(self.stack_std, out_prefix+'_std.tif', ds) iolib.writeGTiff(self.stack_min, out_prefix+'_min.tif', ds) iolib.writeGTiff(self.stack_max, out_prefix+'_max.tif', ds) if self.med: iolib.writeGTiff(self.stack_med, out_prefix+'_med.tif', ds) iolib.writeGTiff(self.stack_nmad, out_prefix+'_nmad.tif', ds)
def write_trend(self): #stat_list = ['stack_trend', 'stack_intercept', 'stack_detrended_std', 'stack_rsquared'] stat_list = ['stack_trend', 'stack_intercept', 'stack_detrended_std'] if any([not hasattr(self, i) for i in stat_list]): self.compute_trend() print("Writing out trend") #Create dummy ds - might want to use vrt here instead driver = gdal.GetDriverByName("MEM") ds = driver.Create('', self.ma_stack.shape[2], self.ma_stack.shape[1], 1, gdal.GDT_Float32) ds.SetGeoTransform(self.gt) ds.SetProjection(self.proj) #Write out with malib, should preserve ma type out_prefix = os.path.splitext(self.stack_fn)[0] iolib.writeGTiff(self.stack_trend, out_prefix+'_trend.tif', ds) iolib.writeGTiff(self.stack_intercept, out_prefix+'_intercept.tif', ds) iolib.writeGTiff(self.stack_detrended_std, out_prefix+'_detrended_std.tif', ds) #iolib.writeGTiff(self.stack_rsquared, out_prefix+'_rsquared.tif', ds)
def mem_ds(res, extent, srs=None, dtype=gdal.GDT_Float32): """Create a new GDAL Dataset in memory Useful for various applications that require a Dataset """ #These round down to int #dst_ns = int((extent[2] - extent[0])/res) #dst_nl = int((extent[3] - extent[1])/res) #This should pad by 1 pixel, but not if extent and res were calculated together to give whole int dst_ns = int((extent[2] - extent[0])/res + 0.99) dst_nl = int((extent[3] - extent[1])/res + 0.99) m_ds = gdal.GetDriverByName('MEM').Create('', dst_ns, dst_nl, 1, dtype) m_gt = [extent[0], res, 0, extent[3], 0, -res] m_ds.SetGeoTransform(m_gt) if srs is not None: m_ds.SetProjection(srs.ExportToWkt()) return m_ds #Modify proj/gt of dst_fn in place
def save_geotiff(self,filename, dtype=gdal.GDT_Float32, **kwargs): """ Save a GeoTIFF of the raster currently loaded. Georeferenced and subset according to the current raster. :params filename: the path and filename to save the file to. :type filename: str :params dtype: GDAL datatype, defaults to Float32. :type dtype: int :returns: 1 on success. :rtype: int """ if self.r is None: raise ValueError('There are no image data loaded. No file can be created.') simple_write_geotiff(filename, self.r, self.trans, wkt=self.srs.ExportToWkt(), dtype=dtype, **kwargs)
def array2Raster(result_arr,affine_par,dstFilename): # save arr to geotiff driver = gdal.GetDriverByName('GTiff') [A, D, B, E, C, F] = affine_par if result_arr.ndim == 2: [height,width] = result_arr.shape dataset = driver.Create(dstFilename, width, height, numBand, gdal.GDT_Float32) dataset.SetGeoTransform((C, A, B, F, D, E)) dataset.GetRasterBand(1).WriteArray(array[:, :]) #dataset.GetRasterBand(1).SetNoDataValue(0) elif result_arr.ndim == 3: [height,width,numBand] = result_arr.shape dataset = driver.Create(dstFilename, width, height, numBand, gdal.GDT_Float32) dataset.SetGeoTransform((C, A, B, F, D, E)) for i in range(numBand): dataset.GetRasterBand(i + 1).WriteArray(result_arr[:, :, i]) #dataset.GetRasterBand(i + 1).SetNoDataValue(-999.0) dataset.FlushCache() # Write to disk
def rasterize_polygons(polygon_shp, template_raster, out_raster, field): """ generate a categorical raster based on polygons :rtype: None :param polygon_shp: input polygon shapefile :param template_raster: raster template for cellsize and extent :param out_raster: output raster file """ # Source: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html gdal.UseExceptions() # Get template raster specs src_ds = gdal.Open(template_raster) if src_ds is None: print 'Unable to open %s' % template_raster sys.exit(1) try: srcband = src_ds.GetRasterBand(1) except RuntimeError, e: print 'No band %i found' % 0 print e sys.exit(1) # Open the data source and read in the extent source_ds = ogr.Open(polygon_shp) source_layer = source_ds.GetLayer() target_ds = gdal.GetDriverByName('GTiff').Create(out_raster, src_ds.RasterXSize, src_ds.RasterYSize, 1, gdal.GDT_Float32) target_ds.SetGeoTransform(src_ds.GetGeoTransform()) target_ds.SetProjection(src_ds.GetProjection()) band = target_ds.GetRasterBand(1) band.SetNoDataValue(srcband.GetNoDataValue()) # Rasterize gdal.RasterizeLayer(target_ds, [1], source_layer, options=["ATTRIBUTE={}".format(field)])
def write_file(output_name,output_array,GeoT,xsize,ysize,proJ,driverName='GTiff'): print "creating", output_name dr=gdal.GetDriverByName(driverName) dr.Register() do=dr.Create(output_name,xsize,ysize,1,gdal.GDT_Float32) do.SetGeoTransform(GeoT) do.SetProjection(proJ) do.GetRasterBand(1).WriteArray(output_array) do.GetRasterBand(1).SetNoDataValue(32767) do=None
def write_file(output_name,output_array,GeoT,xsize,ysize,proJ,driverName='GTiff'): print "creating", output_name dr=gdal.GetDriverByName(driverName) dr.Register() do=dr.Create(output_name,xsize,ysize,1,gdal.GDT_Float32) do.SetGeoTransform(GeoT) do.SetProjection(proJ) do.GetRasterBand(1).WriteArray(output_array) do.GetRasterBand(1).SetNoDataValue(32767) do=None # Function to calculate the moving average
def write_datestack(self): #stat_list = ['dt_stack_ptp', 'dt_stack_mean', 'dt_stack_min', 'dt_stack_max', 'dt_stack_center'] stat_list = ['dt_stack_ptp', 'dt_stack_min', 'dt_stack_max', 'dt_stack_center'] if any([not hasattr(self, i) for i in stat_list]): #self.make_datestack() self.compute_dt_stats() print("Writing out datestack stats") #Create dummy ds - might want to use vrt here instead driver = gdal.GetDriverByName("MEM") ds = driver.Create('', self.dt_stack_ptp.shape[1], self.dt_stack_ptp.shape[0], 1, gdal.GDT_Float32) ds.SetGeoTransform(self.gt) ds.SetProjection(self.proj) #Write out with malib, should preserve ma type out_prefix = os.path.splitext(self.stack_fn)[0] iolib.writeGTiff(self.dt_stack_ptp, out_prefix+'_dt_ptp.tif', ds) self.dt_stack_ptp.set_fill_value(-9999) #iolib.writeGTiff(self.dt_stack_mean, out_prefix+'_dt_mean.tif', ds) #self.dt_stack_mean.set_fill_value(-9999) iolib.writeGTiff(self.dt_stack_min, out_prefix+'_dt_min.tif', ds) self.dt_stack_min.set_fill_value(-9999) iolib.writeGTiff(self.dt_stack_max, out_prefix+'_dt_max.tif', ds) self.dt_stack_max.set_fill_value(-9999) iolib.writeGTiff(self.dt_stack_center, out_prefix+'_dt_center.tif', ds) self.dt_stack_center.set_fill_value(-9999) #Note: want ot change the variable names min/max here #Might be better to save out as multiband GTiff here
def gdaldem_mem_ma(ma, ds=None, res=None, extent=None, srs=None, processing='hillshade', returnma=False): """ Wrapper to allow gdaldem calculations for arbitrary NumPy masked array input Untested, work in progress placeholder Should only need to specify res, can caluclate local gt, cartesian srs """ if ds is None: ds = mem_ds(res, extent, srs=None, dtype=gdal.GDT_Float32) b = ds.GetRasterBand(1) b.WriteArray(ma) return gdaldem_mem_ds(ds, processing=processing, returnma=returnma)
def make_metatile(pile_directory, desired_value): metatile_filename = os.path.join(pile_directory, desired_value + '.tif') metatile_ds = gdal.GetDriverByName('GTiff').Create( metatile_filename, 4096, 4096, 1, gdal.GDT_Float32) # TODO: Try to add georeferencing... return metatile_filename, metatile_ds
def arr2Raster(data_name,Dsm_arr): [h,w] = Dsm_arr.shape driver = gdal.GetDriverByName('GTiff') dataset = driver.Create( data_name, w, h, 1, gdal.GDT_Float32) dataset.GetRasterBand(1).WriteArray(Dsm_arr[:, :]) dataset.FlushCache()
def from_array(cls, raster, geo_transform, proj4, gdal_dtype=gdal.GDT_Float32, nodata=None): """ Create a georaster object from numpy array and georeferencing information. :param raster: 2-D NumPy array of raster to load :type raster: np.array :param geo_transform: a Geographic Transform tuple of the form \ (top left x, w-e cell size, 0, top left y, 0, n-s cell size (-ve)) :type geo_transform: tuple :param proj4: a proj4 string representing the raster projection :type proj4: str :param gdal_dtype: a GDAL data type (default gdal.GDT_Float32) :type gdal_dtype: int :param nodata: None or the nodata value for this array :type nodata: None, int, float, str :returns: GeoRaster object :rtype: GeoRaster """ if len(raster.shape) > 2: nbands = raster.shape[2] else: nbands = 1 # Create a GDAL memory raster to hold the input array mem_drv = gdal.GetDriverByName('MEM') source_ds = mem_drv.Create('', raster.shape[1], raster.shape[0], nbands, gdal_dtype) # Set geo-referencing source_ds.SetGeoTransform(geo_transform) srs = osr.SpatialReference() srs.ImportFromProj4(proj4) source_ds.SetProjection(srs.ExportToWkt()) # Write input array to the GDAL memory raster for b in range(0,nbands): if nbands > 1: r = raster[:,:,b] else: r = raster source_ds.GetRasterBand(b+1).WriteArray(r) if nodata != None: source_ds.GetRasterBand(b+1).SetNoDataValue(nodata) # Return a georaster instance instantiated by the GDAL raster return cls(source_ds)
def reproject_dataset_example(dataset, dataset_example, method=1): # open dataset that must be transformed try: if dataset.split('.')[-1] == 'tif': g = gdal.Open(dataset) else: g = dataset except: g = dataset epsg_from = Get_epsg(g) # open dataset that is used for transforming the dataset try: if dataset_example.split('.')[-1] == 'tif': gland = gdal.Open(dataset_example) else: gland = dataset_example except: gland = dataset_example epsg_to = Get_epsg(gland) # Set the EPSG codes osng = osr.SpatialReference() osng.ImportFromEPSG(epsg_to) wgs84 = osr.SpatialReference() wgs84.ImportFromEPSG(epsg_from) # Get shape and geo transform from example geo_land = gland.GetGeoTransform() col=gland.RasterXSize rows=gland.RasterYSize # Create new raster mem_drv = gdal.GetDriverByName('MEM') dest1 = mem_drv.Create('', col, rows, 1, gdal.GDT_Float32) dest1.SetGeoTransform(geo_land) dest1.SetProjection(osng.ExportToWkt()) # Perform the projection/resampling if method is 1: gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_NearestNeighbour) if method is 2: gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Bilinear) if method is 3: gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Lanczos) if method is 4: gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Average) return(dest1)
def Vector_to_Raster(Dir, shapefile_name, reference_raster_data_name): """ This function creates a raster of a shp file Keyword arguments: Dir -- str: path to the basin folder shapefile_name -- 'C:/....../.shp' str: Path from the shape file reference_raster_data_name -- 'C:/....../.tif' str: Path to an example tiff file (all arrays will be reprojected to this example) """ from osgeo import gdal, ogr geo, proj, size_X, size_Y=Open_array_info(reference_raster_data_name) x_min = geo[0] x_max = geo[0] + size_X * geo[1] y_min = geo[3] + size_Y * geo[5] y_max = geo[3] pixel_size = geo[1] # Filename of the raster Tiff that will be created Dir_Basin_Shape = os.path.join(Dir,'Basin') if not os.path.exists(Dir_Basin_Shape): os.mkdir(Dir_Basin_Shape) Basename = os.path.basename(shapefile_name) Dir_Raster_end = os.path.join(Dir_Basin_Shape, os.path.splitext(Basename)[0]+'.tif') # Open the data source and read in the extent source_ds = ogr.Open(shapefile_name) source_layer = source_ds.GetLayer() # Create the destination data source x_res = int(round((x_max - x_min) / pixel_size)) y_res = int(round((y_max - y_min) / pixel_size)) # Create tiff file target_ds = gdal.GetDriverByName('GTiff').Create(Dir_Raster_end, x_res, y_res, 1, gdal.GDT_Float32, ['COMPRESS=LZW']) target_ds.SetGeoTransform(geo) srse = osr.SpatialReference() srse.SetWellKnownGeogCS(proj) target_ds.SetProjection(srse.ExportToWkt()) band = target_ds.GetRasterBand(1) target_ds.GetRasterBand(1).SetNoDataValue(-9999) band.Fill(-9999) # Rasterize the shape and save it as band in tiff file gdal.RasterizeLayer(target_ds, [1], source_layer, None, None, [1], ['ALL_TOUCHED=TRUE']) target_ds = None # Open array Raster_Basin = Open_tiff_array(Dir_Raster_end) return(Raster_Basin)
def dump_raster(self, filename, driver='GTiff', attr=None, pixel_size=1., remove=True): """ Output layer to GDAL Rasterfile Parameters ---------- filename : string path to shape-filename driver : string GDAL Raster Driver attr : string attribute to burn into raster pixel_size : float pixel Size in source units remove : bool if True removes existing output file """ layer = self.ds.GetLayer() layer.ResetReading() x_min, x_max, y_min, y_max = layer.GetExtent() cols = int((x_max - x_min) / pixel_size) rows = int((y_max - y_min) / pixel_size) # Todo: at the moment, always writing floats ds_out = io.gdal_create_dataset('MEM', '', cols, rows, 1, gdal_type=gdal.GDT_Float32) ds_out.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size)) proj = layer.GetSpatialRef() if proj is None: proj = self._srs ds_out.SetProjection(proj.ExportToWkt()) band = ds_out.GetRasterBand(1) band.FlushCache() print("Rasterize layers") if attr is not None: gdal.RasterizeLayer(ds_out, [1], layer, burn_values=[0], options=["ATTRIBUTE={0}".format(attr), "ALL_TOUCHED=TRUE"], callback=gdal.TermProgress) else: gdal.RasterizeLayer(ds_out, [1], layer, burn_values=[1], options=["ALL_TOUCHED=TRUE"], callback=gdal.TermProgress) io.write_raster_dataset(filename, ds_out, driver, remove=remove) del ds_out
def gdal_create_dataset(drv, name, cols=0, rows=0, bands=0, gdal_type=gdal.GDT_Unknown, remove=False): """Creates GDAL.DataSet object. .. versionadded:: 0.7.0 .. versionchanged:: 0.11.0 - changed parameters to keyword args - added 'bands' as parameter Parameters ---------- drv : string GDAL driver string name : string path to filename cols : int # of columns rows : int # of rows bands : int # of raster bands gdal_type : raster data type eg. gdal.GDT_Float32 remove : bool if True, existing gdal.Dataset will be removed before creation Returns ------- out : gdal.Dataset object """ driver = gdal.GetDriverByName(drv) metadata = driver.GetMetadata() if not metadata.get('DCAP_CREATE', False): raise IOError("Driver %s doesn't support Create() method.".format(drv)) if remove: if os.path.exists(name): driver.Delete(name) ds = driver.Create(name, cols, rows, bands, gdal_type) return ds
def open_data(filename): ''' The function open and load the image given its name. The type of the data is checked from the file and the scipy array is initialized accordingly. Input: filename: the name of the file Output: im: the data cube GeoTransform: the geotransform information Projection: the projection information ''' data = gdal.Open(filename,gdal.GA_ReadOnly) if data is None: print 'Impossible to open '+filename exit() nc = data.RasterXSize nl = data.RasterYSize d = data.RasterCount # Get the type of the data gdal_dt = data.GetRasterBand(1).DataType if gdal_dt == gdal.GDT_Byte: dt = 'uint8' elif gdal_dt == gdal.GDT_Int16: dt = 'int16' elif gdal_dt == gdal.GDT_UInt16: dt = 'uint16' elif gdal_dt == gdal.GDT_Int32: dt = 'int32' elif gdal_dt == gdal.GDT_UInt32: dt = 'uint32' elif gdal_dt == gdal.GDT_Float32: dt = 'float32' elif gdal_dt == gdal.GDT_Float64: dt = 'float64' elif gdal_dt == gdal.GDT_CInt16 or gdal_dt == gdal.GDT_CInt32 or gdal_dt == gdal.GDT_CFloat32 or gdal_dt == gdal.GDT_CFloat64 : dt = 'complex64' else: print 'Data type unkown' exit() # Initialize the array if d == 1: im = sp.empty((nl,nc),dtype=dt) else: im = sp.empty((nl,nc,d),dtype=dt) if d == 1: im[:,:]=data.GetRasterBand(1).ReadAsArray() else : for i in range(d): im[:,:,i]=data.GetRasterBand(i+1).ReadAsArray() GeoTransform = data.GetGeoTransform() Projection = data.GetProjection() data = None return im,GeoTransform,Projection
def write_data(outname,im,GeoTransform,Projection): ''' The function write the image on the hard drive. Input: outname: the name of the file to be written im: the image cube GeoTransform: the geotransform information Projection: the projection information Output: Nothing -- ''' nl = im.shape[0] nc = im.shape[1] if im.ndim == 2: d=1 else: d = im.shape[2] driver = gdal.GetDriverByName('GTiff') dt = im.dtype.name # Get the data type if dt == 'bool' or dt == 'uint8': gdal_dt=gdal.GDT_Byte elif dt == 'int8' or dt == 'int16': gdal_dt=gdal.GDT_Int16 elif dt == 'uint16': gdal_dt=gdal.GDT_UInt16 elif dt == 'int32': gdal_dt=gdal.GDT_Int32 elif dt == 'uint32': gdal_dt=gdal.GDT_UInt32 elif dt == 'int64' or dt == 'uint64' or dt == 'float16' or dt == 'float32': gdal_dt=gdal.GDT_Float32 elif dt == 'float64': gdal_dt=gdal.GDT_Float64 elif dt == 'complex64': gdal_dt=gdal.GDT_CFloat64 else: print 'Data type non-suported' exit() dst_ds = driver.Create(outname,nc,nl, d, gdal_dt) dst_ds.SetGeoTransform(GeoTransform) dst_ds.SetProjection(Projection) if d==1: out = dst_ds.GetRasterBand(1) out.WriteArray(im) out.FlushCache() else: for i in range(d): out = dst_ds.GetRasterBand(i+1) out.WriteArray(im[:,:,i]) out.FlushCache() dst_ds = None
