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)
python类GDT_Float32()的实例源码
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
moving_average_NH.py 文件源码
项目:Global_GPP_VPM_NCEP_C3C4
作者: zhangyaonju
项目源码
文件源码
阅读 19
收藏 0
点赞 0
评论 0
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
moving_average_SH.py 文件源码
项目:Global_GPP_VPM_NCEP_C3C4
作者: zhangyaonju
项目源码
文件源码
阅读 19
收藏 0
点赞 0
评论 0
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_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_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