The following are code examples for showing how to use gdal.GetDriverByName().
Example 1
def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array, nodata, EPSG):
"""This function take a regular array to create a raster with it"""
print("I am dealing with nodata values")
array[np.isnan(array)] = nodata # Dealing with Nodata values
print("I am writing the raster")
cols = array.shape[1]
rows = array.shape[0]
originX = rasterOrigin[0]
originY = rasterOrigin[1]
driver = gdal.GetDriverByName('ENVI')
outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float64)
outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
outband = outRaster.GetRasterBand(1)
#outband.SetNoDataValue(nodata)
outband.WriteArray(array)
#outband.SetNoDataValue(nodata)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromEPSG(EPSG)
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()
Example 2
def writeFile(filename,geotransform,geoprojection,data):
(x,y) = data.shape
format = "GTiff"
noDataValue = -9999
driver = gdal.GetDriverByName(format)
# you can change the dataformat but be sure to be able to store negative values including -9999
dst_datatype = gdal.GDT_Float32
#print(data)
dst_ds = driver.Create(filename,y,x,1,dst_datatype)
dst_ds.GetRasterBand(1).WriteArray(data)
dst_ds.GetRasterBand(1).SetNoDataValue( noDataValue )
dst_ds.SetGeoTransform(geotransform)
dst_ds.SetProjection(geoprojection)
return 1
Example 3
def writeFile(filename,geotransform,geoprojection,data):
(x,y) = data.shape
format = "GTiff"
noDataValue = -9999
driver = gdal.GetDriverByName(format)
# you can change the dataformat but be sure to be able to store negative values including -9999
dst_datatype = gdal.GDT_Float32
#print(data)
dst_ds = driver.Create(filename,y,x,1,dst_datatype)
dst_ds.GetRasterBand(1).WriteArray(data)
dst_ds.GetRasterBand(1).SetNoDataValue( noDataValue )
dst_ds.SetGeoTransform(geotransform)
dst_ds.SetProjection(geoprojection)
return 1
Example 4
def set_parameters(self, tile_params):
self.extent = tile_params.get('extent')
self.tile_width = tile_params.get('width', 256)
self.tile_height = tile_params.get('height', 256)
self.min_zoom = tile_params.get('min_zoom')
self.max_zoom = tile_params.get('max_zoom')
tile_format = tile_params['format']
options = []
if tile_format == 'JPG':
tile_format = 'JPEG'
options = ['QUALITY=%s' % tile_params.get('quality', 75)]
driver = gdal.GetDriverByName('MBTiles')
ds = driver.Create(self.filename, 1, 1, 1, options=['TILE_FORMAT=%s' % tile_format] + options)
ds = None
self._execute_sqlite(
"INSERT INTO metadata(name, value) VALUES ('{}', '{}');".format('minzoom', self.min_zoom),
"INSERT INTO metadata(name, value) VALUES ('{}', '{}');".format('maxzoom', self.max_zoom),
# will be set properly after writing all tiles
"INSERT INTO metadata(name, value) VALUES ('{}', '');".format('bounds')
)
self._zoom = None
Example 5
def write(self, array, name):
self.name = name
driver= gdal.GetDriverByName("GTiff")
if self.DataType == "Float32":
#print gdal.GDT_Float32
self.DataType = gdal.GDT_Float32
self.outdata = array
self.outdata[np.isnan(self.outdata)] = self.NDV
# Create output file
band = 1
DataSet = driver.Create(self.name, self.x, self.y,
band, self.DataType)
DataSet.SetGeoTransform(self.GeoT)
DataSet.SetProjection(self.prj.ExportToWkt())
# Write data array
DataSet.GetRasterBand(1).WriteArray(self.outdata)
DataSet.GetRasterBand(1).SetNoDataValue(self.NDV)
DataSet = None
Example 6
def Get_Extent(input_lyr):
"""
Obtain the input layer extent (xmin, ymin, xmax, ymax)
"""
# Input
filename, ext = os.path.splitext(input_lyr)
if ext.lower() == '.shp':
inp_driver = ogr.GetDriverByName('ESRI Shapefile')
inp_source = inp_driver.Open(input_lyr)
inp_lyr = inp_source.GetLayer()
x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
inp_lyr = None
inp_source = None
elif ext.lower() == '.tif':
inp_lyr = gdal.Open(input_lyr)
inp_transform = inp_lyr.GetGeoTransform()
x_min = inp_transform[0]
x_max = x_min + inp_transform[1] * inp_lyr.RasterXSize
y_max = inp_transform[3]
y_min = y_max + inp_transform[5] * inp_lyr.RasterYSize
inp_lyr = None
else:
raise Exception('The input data type is not recognized')
return (x_min, y_min, x_max, y_max)
Example 7
def save_raster_simple(array, path, dst_filename):
""" Save an array base on an existing raster """
example = gdal.Open(path)
x_pixels = array.shape[1] # number of pixels in x
y_pixels = array.shape[0] # number of pixels in y
bands = 1
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create(dst_filename,x_pixels, y_pixels, bands,
gdal.GDT_Int32)
geotrans=example.GetGeoTransform() #get GeoTranform from existed 'data0'
proj=example.GetProjection() #you can get from a exsited tif or import
dataset.SetGeoTransform(geotrans)
dataset.SetProjection(proj)
dataset.GetRasterBand(1).WriteArray(array[:,:])
dataset.FlushCache()
Example 8
def save_raster(array, path, dst_filename):
""" Save the final multiband array based on an existing raster """
example = gdal.Open(path)
x_pixels = array.shape[2] # number of pixels in x
y_pixels = array.shape[1] # number of pixels in y
bands = array.shape[0]
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create(dst_filename,x_pixels,
y_pixels, bands ,gdal.GDT_Int32)
geotrans=example.GetGeoTransform() #get GeoTranform from existed 'data0'
proj=example.GetProjection() #you can get from a exsited tif or import
dataset.SetGeoTransform(geotrans)
dataset.SetProjection(proj)
for b in range(bands):
dataset.GetRasterBand(b+1).WriteArray(array[b,:,:])
dataset.FlushCache()
Example 9
def save_raster_memory(array, path):
""" Save a raster into memory """
example = gdal.Open(path)
x_pixels = array.shape[1] # number of pixels in x
y_pixels = array.shape[0] # number of pixels in y
driver = gdal.GetDriverByName('MEM')
dataset = driver.Create('',x_pixels, y_pixels, 1,gdal.GDT_Int32)
dataset.GetRasterBand(1).WriteArray(array[:,:])
# follow code is adding GeoTranform and Projection
geotrans=example.GetGeoTransform() #get GeoTranform from existed 'data0'
proj=example.GetProjection() #you can get from a exsited tif or import
dataset.SetGeoTransform(geotrans)
dataset.SetProjection(proj)
return dataset
Example 10
def create_temp_tiff(self, name, content=np.ones([3, 3, 3]), geotransform=(10, 10, 0, 10, 0, -10)):
"""Creates a temporary geotiff in self.path
"""
if len(content.shape) != 3:
raise IndexError
path = os.path.join(self.path, name)
driver = gdal.GetDriverByName('GTiff')
new_image = driver.Create(
path,
xsize=content.shape[1],
ysize=content.shape[2],
bands=content.shape[0],
eType=gdal.GDT_Byte
)
new_image.SetGeoTransform(geotransform)
for band in range(content.shape[0]):
raster_band = new_image.GetRasterBand(band+1)
raster_band.WriteArray(content[band, ...].T)
new_image.SetProjection(self.srs.ExportToWkt())
new_image.FlushCache()
self.images.append(new_image)
self.image_paths.append(path)
Example 11
def create_temp_shape(self, name, point_list):
vector_file = os.path.join(self.temp_dir.name, name)
shape_driver = ogr.GetDriverByName("ESRI Shapefile") # Depreciated; replace at some point
vector_data_source = shape_driver.CreateDataSource(vector_file)
vector_layer = vector_data_source.CreateLayer("geometry", self.srs, geom_type=ogr.wkbPolygon)
ring = ogr.Geometry(ogr.wkbLinearRing)
for point in point_list:
ring.AddPoint(point[0], point[1])
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
vector_feature_definition = vector_layer.GetLayerDefn()
vector_feature = ogr.Feature(vector_feature_definition)
vector_feature.SetGeometry(poly)
vector_layer.CreateFeature(vector_feature)
vector_layer.CreateField(ogr.FieldDefn("class", ogr.OFTInteger))
feature = ogr.Feature(vector_layer.GetLayerDefn())
feature.SetField("class", 3)
vector_data_source.FlushCache()
self.vectors.append(vector_data_source) # Check this is the right thing to be saving here
self.vector_paths.append(vector_file)
Example 12
def geotiff_dir():
"""
Returns
-------
A pointer to a temporary folder that contains a 3-band geotiff
of 3x3, with all values being 1.
"""
tempDir = tempfile.TemporaryDirectory()
fileformat = "GTiff"
driver = gdal.GetDriverByName(fileformat)
metadata = driver.GetMetadata()
tempPath = os.path.join(tempDir.name)
testDataset = driver.Create(os.path.join(tempDir.name, "tempTiff.tif"),
xsize=3, ysize=3, bands=3, eType=gdal.GDT_CFloat32)
for i in range(3):
testDataset.GetRasterBand(i+1).WriteArray(np.ones([3, 3]))
testDataset = None
yield tempPath
tempDir.cleanup()
Example 13
def driver_create(self):
"""
Use gdal.Driver.Create to create out_image.tiff.
"""
file_name = os.path.join(self.data_dir, "out_image.tiff")
image_format = "GTiff"
driver = gdal.GetDriverByName(str(image_format))
data_source = driver.Create(file_name, 50, 50, 1, gdal.GDT_Byte)
raster = np.ones((50, 50), dtype=np.uint8)
raster[10:40, 10:40] = 0
raster = raster * 255
data_source.GetRasterBand(1).WriteArray(raster)
# Avoid PermissionError on Windows when trying to delete
# file_name. From:
# http://stackoverflow.com/questions/22068148/extremely-frustrating-behavior-with-windowserror-error-32-to-remove-temporary
data_source.FlushCache()
driver = None
data_source = None
os.remove(file_name)
Example 14
def save_img(data, geotransform, proj, outPath, noDataValue=np.nan, fieldNames=[]):
# Start the gdal driver for GeoTIFF
if outPath == "MEM":
driver = gdal.GetDriverByName("MEM")
driverOpt = []
shape = data.shape
if len(shape) > 2:
ds = driver.Create(outPath, shape[1], shape[0], shape[2], gdal.GDT_Float32, driverOpt)
ds.SetProjection(proj)
ds.SetGeoTransform(geotransform)
for i in range(shape[2]):
ds.GetRasterBand(i+1).WriteArray(data[:, :, i])
ds.GetRasterBand(i+1).SetNoDataValue(noDataValue)
else:
ds = driver.Create(outPath, shape[1], shape[0], 1, gdal.GDT_Float32, driverOpt)
ds.SetProjection(proj)
ds.SetGeoTransform(geotransform)
ds.GetRasterBand(1).WriteArray(data)
ds.GetRasterBand(1).SetNoDataValue(noDataValue)
return ds
Example 15
def rasterize(data, vectorSrc, field, outFile):
dataSrc = gdal.Open(data)
import ogr
shp = ogr.Open(vectorSrc)
lyr = shp.GetLayer()
driver = gdal.GetDriverByName('GTiff')
dst_ds = driver.Create(
outFile,
dataSrc.RasterXSize,
dataSrc.RasterYSize,
1,
gdal.GDT_UInt16)
dst_ds.SetGeoTransform(dataSrc.GetGeoTransform())
dst_ds.SetProjection(dataSrc.GetProjection())
if field is None:
gdal.RasterizeLayer(dst_ds, [1], lyr, None)
else:
OPTIONS = ['ATTRIBUTE=' + field]
gdal.RasterizeLayer(dst_ds, [1], lyr, None, options=OPTIONS)
data, dst_ds, shp, lyr = None, None, None, None
return outFile
Example 16
def Get_Extent(input_lyr):
"""
Obtain the input layer extent (xmin, ymin, xmax, ymax)
"""
# Input
filename, ext = os.path.splitext(input_lyr)
if ext.lower() == '.shp':
inp_driver = ogr.GetDriverByName('ESRI Shapefile')
inp_source = inp_driver.Open(input_lyr)
inp_lyr = inp_source.GetLayer()
x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
inp_lyr = None
inp_source = None
elif ext.lower() == '.tif':
inp_lyr = gdal.Open(input_lyr)
inp_transform = inp_lyr.GetGeoTransform()
x_min = inp_transform[0]
x_max = x_min + inp_transform[1] * inp_lyr.RasterXSize
y_max = inp_transform[3]
y_min = y_max + inp_transform[5] * inp_lyr.RasterYSize
inp_lyr = None
else:
raise Exception('The input data type is not recognized')
return (x_min, y_min, x_max, y_max)
Example 17
def get_gdal_mem_datset(self): # type: (...) -> gdal.Dataset
driver = gdal.GetDriverByName('MEM')
dtype = self.get_metadata().get_gdal_datatype()
dataset = driver.Create(
'',
self.get_metadata().get_npix_x(),
self.get_metadata().get_npix_y(),
self.get_metadata().get_n_bands(),
dtype)
dataset.SetGeoTransform(self.get_point_calculator().get_geot())
dataset.SetProjection(self.get_point_calculator().get_gdal_projection_wkt())
if self.get_image_data() is None:
self.set_image_data(self.read_all_image_data_from_disk())
for band in range(self.get_metadata().get_n_bands()):
dataset.GetRasterBand(band + 1).WriteArray(self.get_image_band(band))
if self.get_metadata().get_nodata_val() is not None:
dataset.GetRasterBand(band + 1).SetNoDataValue(int(self.get_metadata().get_nodata_val()))
return dataset
Example 18
def write_to_disk(self,
fname, # type: str
compression="LZW" # type: str
): # type: (...) -> None
driver = gdal.GetDriverByName('GTiff')
dtype = self.get_metadata().get_gdal_datatype()
ops = ['COMPRESS=' + compression, "INTERLEAVE=BAND", "TILED=YES"]
dataset = driver.Create(
fname,
self.get_metadata().get_npix_x(),
self.get_metadata().get_npix_y(),
self.get_metadata().get_n_bands(),
dtype,
ops)
dataset.SetGeoTransform(self.get_point_calculator().get_geot())
dataset.SetProjection(self.get_point_calculator().get_gdal_projection_wkt())
if self.get_image_data() is None:
self.set_image_data(self.read_all_image_data_from_disk())
for band in range(self.get_metadata().get_n_bands()):
dataset.GetRasterBand(band + 1).WriteArray(self.get_image_band(band))
if self.get_metadata().get_nodata_val() is not None:
dataset.GetRasterBand(band + 1).SetNoDataValue(int(self.get_metadata().get_nodata_val()))
dataset.FlushCache()
dataset = None
Example 19
def convertToRaster(self,inFName,outFName,xRecs,yRecs,transformArray,keyField=None):
if self.debug == True:
QgsMessageLog.logMessage(self.myself())
# Open the data source and read in the extent
sDs = ogr.Open(inFName)
sLyr = sDs.GetLayer()
gtDriver = gdal.GetDriverByName('GTiff')
tDs = gtDriver.Create(outFName, int(xRecs), int(yRecs), 1, gdal.GDT_Int32)
tDs.SetGeoTransform(transformArray)
# set options
if keyField == None:
#options = ['ALL_TOUCHED=True']
options = []
else:
#options = ['ATTRIBUTE=%s' % keyField,'ALL_TOUCHED=TRUE']
options = ['ATTRIBUTE=%s' % keyField]
# rasterize
gdal.RasterizeLayer(tDs, [1], sLyr, None, None, [1], options)
sLyr = None
sDs = None
tDs = None
Example 20
def create_temp_tiff(self, name, content=np.ones([3, 3, 3]), geotransform=(100, 10, 0, 100, 0, 10)):
"""Creates a temporary geotiff in self.path
"""
if len(content.shape) != 3:
raise IndexError
path = os.path.join(self.path, name)
driver = gdal.GetDriverByName('GTiff')
new_image = driver.Create(
path,
xsize=content.shape[1],
ysize=content.shape[2],
bands=content.shape[0],
eType=gdal.GDT_Byte
)
new_image.SetGeoTransform(geotransform)
for band in range(content.shape[0]):
raster_band = new_image.GetRasterBand(band+1)
raster_band.WriteArray(content[band, ...].T)
new_image.SetProjection(self.srs.ExportToWkt())
new_image.FlushCache()
self.images.append(new_image)
self.image_paths.append(path)
Example 21
def create_temp_shp(self, name):
vector_file = os.path.join(self.temp_dir.name, name)
shape_driver = ogr.GetDriverByName("ESRI Shapefile") # Depreciated; replace at some point
vector_data_source = shape_driver.CreateDataSource(vector_file)
vector_layer = vector_data_source.CreateLayer("geometry", self.srs, geom_type=ogr.wkbPolygon)
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(100.0, 100.0)
ring.AddPoint(100.0, 110.0)
ring.AddPoint(110.0, 110.0)
ring.AddPoint(110.0, 100.0)
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
vector_feature_definition = vector_layer.GetLayerDefn()
vector_feature = ogr.Feature(vector_feature_definition)
vector_feature.SetGeometry(poly)
vector_layer.CreateFeature(vector_feature)
vector_data_source.FlushCache()
self.vectors.append(vector_data_source) # Check this is the right thing to be saving here
self.vector_paths.append(vector_file)
Example 22
def geotiff_dir():
"""
Returns
-------
A pointer to a temporary folder that contains a 3-band geotiff
of 3x3, with all values being 1.
"""
tempDir = tempfile.TemporaryDirectory()
fileformat = "GTiff"
driver = gdal.GetDriverByName(fileformat)
metadata = driver.GetMetadata()
tempPath = os.path.join(tempDir.name)
testDataset = driver.Create(os.path.join(tempDir.name, "tempTiff.tif"),
xsize=3, ysize=3, bands=3, eType=gdal.GDT_CFloat32)
for i in range(3):
testDataset.GetRasterBand(i+1).WriteArray(np.ones([3, 3]))
testDataset = None
yield tempPath
tempDir.cleanup()
Example 23
def write(self, array, name):
self.name = name
driver= gdal.GetDriverByName("GTiff")
if self.DataType == "Float32":
#print gdal.GDT_Float32
self.DataType = gdal.GDT_Float32
self.outdata = array
self.outdata[np.isnan(self.outdata)] = self.NDV
# Create output file
band = 1
DataSet = driver.Create(self.name, self.x, self.y,
band, self.DataType)
DataSet.SetGeoTransform(self.GeoT)
DataSet.SetProjection(self.prj.ExportToWkt())
# Write data array
DataSet.GetRasterBand(1).WriteArray(self.outdata)
DataSet.GetRasterBand(1).SetNoDataValue(self.NDV)
DataSet = None
Example 24
def SaveRaster(raster,path):
"""
===================================================================
SaveRaster(raster,path)
===================================================================
this function saves a raster to a path
inputs:
----------
1- raster:
[gdal object]
2- path:
[string] a path includng the name of the raster and extention like
path="data/cropped.tif"
Outputs:
----------
the function does not return and data but only save the raster to the hard drive
EX:
----------
SaveRaster(raster,output_path)
"""
#### input data validation
# data type
assert type(raster)==gdal.Dataset, "src should be read using gdal (gdal dataset please read it using gdal library) "
assert type(path)== str, "Raster_path input should be string type"
# input values
ext=path[-4:]
assert ext == ".tif", "please add the extension at the end of the path input"
driver = gdal.GetDriverByName ( "GTiff" )
dst_ds = driver.CreateCopy( path, raster, 0 )
dst_ds = None # Flush the dataset to disk
Example 25
def SaveRaster(raster,path):
"""
===================================================================
SaveRaster(raster,path)
===================================================================
this function saves a raster to a path
inputs:
----------
1- raster:
[gdal object]
2- path:
[string] a path includng the name of the raster and extention like
path="data/cropped.tif"
Outputs:
----------
the function does not return and data but only save the raster to the hard drive
EX:
----------
SaveRaster(raster,output_path)
"""
#### input data validation
# data type
assert type(raster)==gdal.Dataset, "src should be read using gdal (gdal dataset please read it using gdal library) "
assert type(path)== str, "Raster_path input should be string type"
# input values
ext=path[-4:]
assert ext == ".tif", "please add the extension at the end of the path input"
driver = gdal.GetDriverByName ( "GTiff" )
dst_ds = driver.CreateCopy( path, raster, 0 )
dst_ds = None # Flush the dataset to disk
Example 26
def Open_bil_array(bil_filename, band = 1):
"""
Opening a bil array.
Keyword Arguments:
bil_filename -- 'C:/file/to/path/file.bil'
string that defines the input tiff file or gdal file
band -- integer
Defines the band of the tiff that must be opened.
"""
gdal.GetDriverByName('EHdr').Register()
img = gdal.Open(bil_filename)
Data = img.GetRasterBand(band).ReadAsArray()
return(Data)
Example 27
def List_Fields(input_lyr):
"""
Lists the field names of input layer
"""
# Input
if isinstance(input_lyr, str):
inp_driver = ogr.GetDriverByName('ESRI Shapefile')
inp_source = inp_driver.Open(input_lyr, 0)
inp_lyr = inp_source.GetLayer()
inp_lyr_defn = inp_lyr.GetLayerDefn()
elif isinstance(input_lyr, ogr.Layer):
inp_lyr_defn = input_lyr.GetLayerDefn()
# List
names_ls = []
# Loop
for j in range(0, inp_lyr_defn.GetFieldCount()):
field_defn = inp_lyr_defn.GetFieldDefn(j)
names_ls.append(field_defn.GetName())
# Save and/or close the data sources
inp_source = None
# Return
return names_ls
Example 28
def Array_to_Raster(input_array, output_tiff, ll_corner, cellsize,
srs_wkt):
"""
Saves an array into a raster file
"""
# Output
out_driver = gdal.GetDriverByName('GTiff')
if os.path.exists(output_tiff):
out_driver.Delete(output_tiff)
y_ncells, x_ncells = input_array.shape
gdal_datatype = gdaltype_from_dtype(input_array.dtype)
out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
1, gdal_datatype)
out_band = out_source.GetRasterBand(1)
out_band.SetNoDataValue(-9999)
out_top_left_x = ll_corner[0]
out_top_left_y = ll_corner[1] + cellsize*y_ncells
out_source.SetGeoTransform((out_top_left_x, cellsize, 0,
out_top_left_y, 0, -cellsize))
out_source.SetProjection(str(srs_wkt))
out_band.WriteArray(input_array)
# Save and/or close the data sources
out_source = None
# Return
return output_tiff
Example 29
def Extract_Band(input_tiff, output_tiff, band_number=1):
"""
Extract and save a raster band into a new raster
"""
# Input
inp_lyr = gdal.Open(input_tiff)
inp_srs = inp_lyr.GetProjection()
inp_transform = inp_lyr.GetGeoTransform()
inp_band = inp_lyr.GetRasterBand(band_number)
inp_array = inp_band.ReadAsArray()
inp_data_type = inp_band.DataType
NoData_value = inp_band.GetNoDataValue()
x_ncells = inp_lyr.RasterXSize
y_ncells = inp_lyr.RasterYSize
# Output
out_driver = gdal.GetDriverByName('GTiff')
if os.path.exists(output_tiff):
out_driver.Delete(output_tiff)
out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
1, inp_data_type)
out_band = out_source.GetRasterBand(1)
out_band.SetNoDataValue(NoData_value)
out_source.SetGeoTransform(inp_transform)
out_source.SetProjection(inp_srs)
out_band.WriteArray(inp_array)
# Save and/or close the data sources
inp_lyr = None
out_source = None
# Return
return output_tiff
Example 30
def do_proximity(config, image, srcarray, label, _class):
""" Get proximity of each pixel to to certain value """
#First create a band in memory that's that's just 1s and 0s
src_ds = gdal.Open( image, gdal.GA_ReadOnly )
srcband = src_ds.GetRasterBand(1)
mem_rast = save_raster_memory(srcarray, image)
mem_band = mem_rast.GetRasterBand(1)
#Now the code behind gdal_sieve.py
maskband = None
drv = gdal.GetDriverByName('GTiff')
#Need to output sieved file
dst_path = config['postprocessing']['deg_class']['prox_dir']
dst_filename = dst_path + '/' + image.split('/')[-1].split('.')[0] + '_' + _class + '.tif'
dst_ds = drv.Create( dst_filename,src_ds.RasterXSize, src_ds.RasterYSize,1,
srcband.DataType )
wkt = src_ds.GetProjection()
if wkt != '':
dst_ds.SetProjection( wkt )
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dstband = dst_ds.GetRasterBand(1)
# Parameters
prog_func = None
options = []
options.append( 'VALUES=' + str(label) )
result = gdal.ComputeProximity(mem_band, dstband, options)#,
#callback = prog_func )
proximity = dstband.ReadAsArray()
return proximity
Example 31
def create_matching_dataset(in_dataset, out_path,
format="GTiff", bands=1, datatype = None):
"""
Creates an empty gdal dataset with the same dimensions, projection and geotransform as in_dataset.
Defaults to 1 band.
Datatype is set from the first layer of in_dataset if unspecified
Parameters
----------
in_dataset
A gdal.Dataset object
out_path
The path to save the copied dataset to
format
The Ggal image format. Defaults to geotiff ("GTiff"); for a full list, see https://gdal.org/drivers/raster/index.html
bands
The number of bands in the dataset. Defaults to one.
datatype
The datatype of the returned dataset. See the introduction for this module.
Returns
-------
An gdal.Dataset of the new, empty dataset that is ready for writing.
"""
driver = gdal.GetDriverByName(format)
if datatype is None:
datatype = in_dataset.GetRasterBand(1).DataType
out_dataset = driver.Create(out_path,
xsize=in_dataset.RasterXSize,
ysize=in_dataset.RasterYSize,
bands=bands,
eType=datatype)
out_dataset.SetGeoTransform(in_dataset.GetGeoTransform())
out_dataset.SetProjection(in_dataset.GetProjection())
return out_dataset
Example 32
def save_array_as_image(array, path, geotransform, projection, format = "GTiff"):
"""
Saves a given array as a geospatial image to disk in the format 'format'. The datatype will be of one corresponding
to
Array must be gdal format: [bands, y, x].
Parameters
----------
array
A Numpy array containing the values to be saved to a raster
path
The path to the location to save the output raster to
geotransform
The geotransform of the image to be saved. See note.
projection
The projection, as wkt, of the image to be saved. See note.
format
Returns
-------
"""
driver = gdal.GetDriverByName(format)
type_code = gdal_array.NumericTypeCodeToGDALTypeCode(array.dtype)
out_dataset = driver.Create(
path,
xsize=array.shape[2],
ysize=array.shape[1],
bands=array.shape[0],
eType=type_code
)
out_dataset.SetGeoTransform(geotransform)
out_dataset.SetProjection(projection)
out_array = out_dataset.GetVirtualMemArray(eAccess=gdal.GA_Update)
out_array[...] = array
out_array = None
out_dataset = None
return path
Example 33
def create_100x100_shp(self, name):
"""Cretes a shapefile with a vector layer named "geometry" containing a 100mx100m square , top left corner
being at wgs coords 10,10.
This polygon has a field, 'class' with a value of 3. Left in for back-compatability"""
# TODO Generalise this
vector_file = os.path.join(self.temp_dir.name, name)
shape_driver = ogr.GetDriverByName("ESRI Shapefile") # Depreciated; replace at some point
vector_data_source = shape_driver.CreateDataSource(vector_file)
vector_layer = vector_data_source.CreateLayer("geometry", self.srs, geom_type=ogr.wkbPolygon)
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(10.0, 10.0)
ring.AddPoint(10.0, 110.0)
ring.AddPoint(110.0, 110.0)
ring.AddPoint(110.0, 10.0)
ring.AddPoint(10.0, 10.0)
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
vector_feature_definition = vector_layer.GetLayerDefn()
vector_feature = ogr.Feature(vector_feature_definition)
vector_feature.SetGeometry(poly)
vector_layer.CreateFeature(vector_feature)
vector_layer.CreateField(ogr.FieldDefn("class", ogr.OFTInteger))
feature = ogr.Feature(vector_layer.GetLayerDefn())
feature.SetField("class", 3)
vector_data_source.FlushCache()
self.vectors.append(vector_data_source) # Check this is the right thing to be saving here
self.vector_paths.append(vector_file)
Example 34
def save_tif(self, mask, file_path):
'''Save the given mask as a .tif file.
Parameters
mask - A mask generated with masker.
file_path - Path of .tif file.
'''
import gdal
driver = gdal.GetDriverByName('GTiff')
x_pixels = mask.shape[1]
y_pixels = mask.shape[0]
dataset = driver.Create(file_path, x_pixels, y_pixels, 1, gdal.GDT_Int32)
if self.file_path is not None:
extension = os.path.splitext(self.file_path)[1].lower()
if extension == '.hdf':
hdfdataset = gdal.Open(self.file_path)
subdataset = hdfdataset.GetSubDatasets()[0][0]
bandfile = gdal.Open(subdataset)
else:
bandfile = gdal.Open(self.file_path)
dataset.SetGeoTransform(bandfile.GetGeoTransform())
dataset.SetProjection(bandfile.GetProjectionRef())
dataset.GetRasterBand(1).WriteArray(mask)
dataset.FlushCache()
Example 35
def __call__(self, rasterBob, spatialReference, filePath = None):
#File will be saved to current working directory if no file path is defined
if filePath == None:
filePath = str(os.getcwd())+"/TestResults"
#Adds the name of the file to be created to the file path
filePath = filePath + "/Results" + str(time.strftime("%m")) + "-" + str(time.strftime("%d"))+".tif"
driver = gdal.GetDriverByName('GTiff') #Obtains the GeoTIFF Driver
#Creates a new .tif file at the indicated file path, with an xSize of w and a ySize of h
newGTiff = driver.Create(filePath, rasterBob.ncols, rasterBob.nrows, 1, gdal.GDT_Float32)
band = newGTiff.GetRasterBand(1)
band.WriteArray(rasterBob.data) #Writes the data from the raster to the new .tif file
band.SetNoDataValue(-1) #Can be changed to whatever value fits the circumstance (still need to decide on the default value)
#Sets up the GeoTransform and the Projection for the file
xRes = rasterBob.w/float(rasterBob.ncols)
yRes = rasterBob.h/float(rasterBob.nrows)
geoTransform = (rasterBob.x, xRes, 0, rasterBob.y, 0, yRes)
newGTiff.SetGeoTransform(geoTransform)
newGTiff.SetProjection(spatialReference.ExportToWkt())
band.FlushCache() #'Saves' the data written to the file
#Closes out the file and raster band
band = None
newGTiff = None
return
Example 36
def __call__(self, STCube, spatialRef, filePath = None):
for time in STCube.timelist:
#File will be saved to current working directory if no file path is defined
if filePath == None:
filePath = str(os.getcwd())+"/TestResults-"+str(time)
#Adds the name of the file to be created to the file path
filePath = filePath + "/Results" + str(time.strftime("%m")) + "-" + str(time.strftime("%d"))+"-"+str(time)+".tif"
driver = gdal.GetDriverByName('GTiff') #Obtains the GeoTIFF Driver
#Creates a new .tif file at the indicated file path, with an xSize of w and a ySize of h
newGTiff = driver.Create(filePath, rasterBob.ncols, rasterBob.nrows, 1, gdal.GDT_Float32)
band = newGTiff.GetRasterBand(1)
band.WriteArray(rasterBob.data) #Writes the data from the raster to the new .tif file
band.SetNoDataValue(-1) #Can be changed to whatever value fits the circumstance (still need to decide on the default value)
#Sets up the GeoTransform and the Projection for the file
xRes = rasterBob.w/float(rasterBob.ncols)
yRes = rasterBob.h/float(rasterBob.nrows)
geoTransform = (rasterBob.x, xRes, 0, rasterBob.y, 0, yRes)
newGTiff.SetGeoTransform(geoTransform)
newGTiff.SetProjection(spatialReference.ExportToWkt())
band.FlushCache() #'Saves' the data written to the file
#Closes out the file and raster band
band = None
newGTiff = None
return
Example 37
def driver_createcopy(self):
"""
Use gdal.Open to load image.tiff and gdal.Driver.CreateCopy to
copy this to out_image.tiff.
"""
file_name = os.path.join(self.data_dir, "image.tiff")
out_file_name = os.path.join(self.data_dir, "out_image.tiff")
data_source = gdal.Open(file_name)
image_format = "GTiff"
driver = gdal.GetDriverByName(str(image_format))
driver.CreateCopy(out_file_name, data_source, 0)
os.remove(out_file_name)
Example 38
def save_to_gtiff(self, fname, left, right):
"""
Save the result image to a geotiff.
Parameters
----------
fname: str
Filename for the anaglyph image.
left: array
Left eye values.
right: array
Right eye values.
"""
driver = gdal.GetDriverByName('GTiff')
ds = driver.Create(fname, self.img_ds.RasterXSize,
self.img_ds.RasterYSize, 3, self.img_rb.DataType)
ds.SetProjection(self.img_ds.GetProjection())
ds.SetGeoTransform(self.img_ds.GetGeoTransform())
# fill out the individual bands
ds.GetRasterBand(1).SetNoDataValue(self.nodata_val)
ds.GetRasterBand(2).SetNoDataValue(self.nodata_val)
ds.GetRasterBand(3).SetNoDataValue(self.nodata_val)
ds.GetRasterBand(1).WriteArray(left)
ds.GetRasterBand(2).WriteArray(right)
ds.GetRasterBand(3).WriteArray(right)
# destroy the dataset object
ds = None
stdout.write(fname + '\n')
Example 39
def saveRaster(nomSortie, image, GeoTransform, Projection):
"""
The function saves the resultant image in the hard disk. Adaptation de write_data sous rasterTool.py
INPUT image : tableau recalculé à 6 bandes.
GeoTransform,Projection : informations from original image
"""
nl, nc, d = np.shape(image)
# ou image.shape; ?
# we create an empty image in GeoTiff
driver = gdal.GetDriverByName('GTiff')
dt = image.dtype.name
if(dt == 'float64'):
gdal_dt = gdal.GDT_Float64
else:
print("Erreur de type de données")
print(dt)
exit()
# blehbleh . "on récupére l'image dans le tableau outname avec ses différents bandes et
# les métadonnées et la projection données en paramétre de la fonction"
dst_ds = driver.Create(nomSortie, nc, nl, d, gdal_dt)
dst_ds.SetGeoTransform(GeoTransform)
dst_ds.SetProjection(Projection)
for i in range(d):
out = dst_ds.GetRasterBand(i + 1)
# on affecte l'image dans le fichier enregistré
out.WriteArray(image[:, :, i])
out.FlushCache()
dst_ds = None
Example 40
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
gdal_dt = getGDALGDT(dt)
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
Example 41
def create_empty_tiff(outname, im, d, GeoTransform, Projection):
'''!@brief Write an empty 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]
driver = gdal.GetDriverByName('GTiff')
dt = im.dtype.name
# Get the data type
gdal_dt = getGDALGDT(dt)
dst_ds = driver.Create(outname, nc, nl, d, gdal_dt)
dst_ds.SetGeoTransform(GeoTransform)
dst_ds.SetProjection(Projection)
return dst_ds
'''
Old function that cannot manage to write on each band outside the script
'''
# 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
Example 42
def create_uniquevalue_tiff(
outname, im, d, GeoTransform, Projection, wholeValue=1, gdal_dt=False):
'''!@brief Write an empty 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]
driver = gdal.GetDriverByName('GTiff')
# Get the data type
if not gdal_dt:
gdal_dt = gdal.GDT_Byte
dst_ds = driver.Create(outname, nc, nl, d, gdal_dt)
dst_ds.SetGeoTransform(GeoTransform)
dst_ds.SetProjection(Projection)
if d == 1:
im[:] = wholeValue
out = dst_ds.GetRasterBand(1)
out.WriteArray(im)
out.FlushCache()
else:
for i in range(d):
im[:, :, i] = wholeValue
out = dst_ds.GetRasterBand(i + 1)
out.WriteArray(im[:, :, i])
out.FlushCache()
dst_ds = None
return outname
Example 43
def List_Fields(input_lyr):
"""
Lists the field names of input layer
"""
# Input
if isinstance(input_lyr, str):
inp_driver = ogr.GetDriverByName('ESRI Shapefile')
inp_source = inp_driver.Open(input_lyr, 0)
inp_lyr = inp_source.GetLayer()
inp_lyr_defn = inp_lyr.GetLayerDefn()
elif isinstance(input_lyr, ogr.Layer):
inp_lyr_defn = input_lyr.GetLayerDefn()
# List
names_ls = []
# Loop
for j in range(0, inp_lyr_defn.GetFieldCount()):
field_defn = inp_lyr_defn.GetFieldDefn(j)
names_ls.append(field_defn.GetName())
# Save and/or close the data sources
inp_source = None
# Return
return names_ls
Example 44
def Array_to_Raster(input_array, output_tiff, ll_corner, cellsize,
srs_wkt):
"""
Saves an array into a raster file
"""
# Output
out_driver = gdal.GetDriverByName('GTiff')
if os.path.exists(output_tiff):
out_driver.Delete(output_tiff)
y_ncells, x_ncells = input_array.shape
gdal_datatype = gdaltype_from_dtype(input_array.dtype)
out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
1, gdal_datatype)
out_band = out_source.GetRasterBand(1)
out_band.SetNoDataValue(-9999)
out_top_left_x = ll_corner[0]
out_top_left_y = ll_corner[1] + cellsize*y_ncells
out_source.SetGeoTransform((out_top_left_x, cellsize, 0,
out_top_left_y, 0, -cellsize))
out_source.SetProjection(str(srs_wkt))
out_band.WriteArray(input_array)
# Save and/or close the data sources
out_source = None
# Return
return output_tiff
Example 45
def Extract_Band(input_tiff, output_tiff, band_number=1):
"""
Extract and save a raster band into a new raster
"""
# Input
inp_lyr = gdal.Open(input_tiff)
inp_srs = inp_lyr.GetProjection()
inp_transform = inp_lyr.GetGeoTransform()
inp_band = inp_lyr.GetRasterBand(band_number)
inp_array = inp_band.ReadAsArray()
inp_data_type = inp_band.DataType
NoData_value = inp_band.GetNoDataValue()
x_ncells = inp_lyr.RasterXSize
y_ncells = inp_lyr.RasterYSize
# Output
out_driver = gdal.GetDriverByName('GTiff')
if os.path.exists(output_tiff):
out_driver.Delete(output_tiff)
out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
1, inp_data_type)
out_band = out_source.GetRasterBand(1)
out_band.SetNoDataValue(NoData_value)
out_source.SetGeoTransform(inp_transform)
out_source.SetProjection(inp_srs)
out_band.WriteArray(inp_array)
# Save and/or close the data sources
inp_lyr = None
out_source = None
# Return
return output_tiff
Example 46
def dump_data(self, timestep, x_analysis, P_analysis, P_analysis_inv,
state_mask, n_params):
drv = gdal.GetDriverByName(self.fmt)
for ii, param in enumerate(self.parameter_list):
if self.prefix is None:
fname = os.path.join(self.folder, "%s_%s.tif" %
(param, timestep.strftime("A%Y%j")))
else:
fname = os.path.join(self.folder, "%s_%s_%s.tif" %
(param, timestep.strftime("A%Y%j"),
self.prefix))
dst_ds = drv.Create(fname, state_mask.shape[1],
state_mask.shape[0], 1,
gdal.GDT_Float32, ['COMPRESS=DEFLATE',
'BIGTIFF=YES',
'PREDICTOR=1',
'TILED=YES'])
dst_ds.SetProjection(self.projection)
dst_ds.SetGeoTransform(self.geotransform)
A = np.zeros(state_mask.shape, dtype=np.float32)
A[state_mask] = x_analysis[ii::n_params]
dst_ds.GetRasterBand(1).WriteArray(A)
for ii, param in enumerate(self.parameter_list):
if self.prefix is None:
fname = os.path.join(self.folder, "%s_%s_unc.tif" %
(param, timestep.strftime("A%Y%j")))
else:
fname = os.path.join(self.folder, "%s_%s_%s_unc.tif" %
(param, timestep.strftime("A%Y%j"),
self.prefix))
dst_ds = drv.Create(fname, state_mask.shape[1],
state_mask.shape[0], 1,
gdal.GDT_Float32, ['COMPRESS=DEFLATE',
'BIGTIFF=YES',
'PREDICTOR=1', 'TILED=YES'])
dst_ds.SetProjection(self.projection)
dst_ds.SetGeoTransform(self.geotransform)
A = np.zeros(state_mask.shape, dtype=np.float32)
A[state_mask] = 1./np.sqrt(P_analysis_inv.diagonal()[ii::n_params])
dst_ds.GetRasterBand(1).WriteArray(A)
Example 47
def vectorToRaster(self,inFName,outFName,xRecs,yRecs,transformArray,keyField=None):
if self.debug == True:
QgsMessageLog.logMessage(self.myself())
# Open the data source and read in the extent
sDs = ogr.Open(inFName)
sLyr = sDs.GetLayer()
gtDriver = gdal.GetDriverByName('GTiff')
tDs = gtDriver.Create(outFName, int(xRecs), int(yRecs), 1, gdal.GDT_Int32)
tDs.SetGeoTransform(transformArray)
# set options
if keyField == None:
#options = ['ALL_TOUCHED=True']
options = []
else:
#options = ['ATTRIBUTE=%s' % keyField,'ALL_TOUCHED=TRUE']
options = ['ATTRIBUTE=%s' % keyField]
# rasterize
#QgsMessageLog.logMessage('rasterize')
gdal.RasterizeLayer(tDs, [1], sLyr, None, None, [1], options)
sLyr = None
sDs = None
tDs = None
#
# Measuring of Raw Values
#
#
# extract records from pu field
Example 48
def SaveGeoCapture(self, task, image, output, p, geotransform):
''' Save Current GeoReferenced Frame '''
ext = ".tiff"
t = "out_" + p + ext
name = "g_" + p
src_file = os.path.join(output, t)
image.save(src_file)
# Opens source dataset
src_ds = gdal.OpenEx(src_file, gdal.OF_RASTER |
gdal.OF_READONLY, open_options=['NUM_THREADS=ALL_CPUS'])
# Open destination dataset
dst_filename = os.path.join(output, name + ext)
dst_ds = gdal.GetDriverByName("GTiff").CreateCopy(dst_filename, src_ds, 0,
options=['TILED=NO', 'BIGTIFF=NO', 'COMPRESS_OVERVIEW=DEFLATE', 'COMPRESS=LZW', 'NUM_THREADS=ALL_CPUS', 'predictor=2'])
src_ds = None
# Get raster projection
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
# Set projection
dst_ds.SetProjection(srs.ExportToWkt())
# Set location
dst_ds.SetGeoTransform(geotransform)
dst_ds.GetRasterBand(1).SetNoDataValue(0)
dst_ds.FlushCache()
# Close files
dst_ds = None
os.remove(src_file)
if task.isCanceled():
return None
return {'task': task.description(),
'file': dst_filename}
Example 49
def _dbf2DF(dbfile, upper=True): #Reads in DBF files and returns Pandas DF
db = io.open(dbfile) #Pysal to open DBF
d = {col: db.by_col(col) for col in db.header} #Convert dbf to dictionary
#pandasDF = pd.DataFrame(db[:]) #Convert to Pandas DF
pandasDF = pd.DataFrame(d) #Convert to Pandas DF
if upper == True: #Make columns uppercase if wanted
pandasDF.columns = map(str.upper, db.header)
db.close()
return pandasDF
##### make a new vector to be written for reference
# outShapefile = outShp
# outDriver = ogr.GetDriverByName("ESRI Shapefile")
#
# # Remove output shapefile if it already exists
# if os.path.exists(outShapefile):
# outDriver.DeleteDataSource(outShapefile)
#
# # Create the output shapefile
# outDataSource = outDriver.CreateDataSource(outShapefile)
# outLayer = outDataSource.CreateLayer("OutLyr", geom_type=ogr.wkbMultiPolygon)
#
# # Add an ID field
# idField = ogr.FieldDefn("id", ogr.OFTInteger)
# outLayer.CreateField(idField)
Example 50
def rasterize_shapefile_like(shpfile, model_raster_fname, dtype, options,
nodata_val=0,):
"""
Given a shapefile, rasterizes it so it has
the exact same extent as the given model_raster
`dtype` is a gdal type like gdal.GDT_Byte
`options` should be a list that will be passed to GDALRasterizeLayers
papszOptions, like ["ATTRIBUTE=vegetation","ALL_TOUCHED=TRUE"]
"""
model_dataset = gdal.Open(model_raster_fname)
shape_dataset = ogr.Open(shpfile)
shape_layer = shape_dataset.GetLayer()
mem_drv = gdal.GetDriverByName('MEM')
mem_raster = mem_drv.Create(
'',
model_dataset.RasterXSize,
model_dataset.RasterYSize,
1,
dtype
)
mem_raster.SetProjection(model_dataset.GetProjection())
mem_raster.SetGeoTransform(model_dataset.GetGeoTransform())
mem_band = mem_raster.GetRasterBand(1)
mem_band.Fill(nodata_val)
mem_band.SetNoDataValue(nodata_val)
# http://gdal.org/gdal__alg_8h.html#adfe5e5d287d6c184aab03acbfa567cb1
# http://gis.stackexchange.com/questions/31568/gdal-rasterizelayer-doesnt-burn-all-polygons-to-raster
err = gdal.RasterizeLayer(
mem_raster,
[1],
shape_layer,
None,
None,
[1],
options
)
assert err == gdal.CE_None
return mem_raster.ReadAsArray()
1692
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
