def save(self, filename, driver=None, *args, **kwargs):
if driver is None:
ext = os.path.splitext(filename)[1][1:]
if ext.lower() in ['tif', 'tiff']:
driver = gdal.GetDriverByName('GTiff')
else:
raise Exception('Unkown extension. Can not determine driver')
bands = self.array.shape[2] if len(self.array.shape)>2 else 1
self.object = driver.Create(filename, self.array.shape[1],
self.array.shape[0], bands, GdalWriter.gdal_array_types[np.dtype(self.dtype)])
if bands==1:
self.object.GetRasterBand(1).WriteArray(self.array)
else:
for band in range(bands):
self.object.GetRasterBand(band+1).WriteArray(self.array[:,:,band])
#del self.object
#Need to be deleted to actually save
# dst_ds.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] )
# srs = osr.SpatialReference()
# srs.SetUTM( 11, 1 )
# srs.SetWellKnownGeogCS( 'NAD27' )
# dst_ds.SetProjection( srs.ExportToWkt() )
# @@@ Common feel functions @@@
python类GetDriverByName()的实例源码
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 make_metatile(pile_directory):
mosaic_filename = os.path.join(pile_directory,'mosaic.tif')
mosaic_ds = gdal.GetDriverByName('GTiff').Create(
mosaic_filename, 4096, 4096, 1, gdal.GDT_UInt16)
# TODO: Try to add georeferencing...
return mosaic_filename, mosaic_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 open_vector(filename, driver=None):
"""Open vector file, return gdal.Dataset and OGR.Layer
.. warning:: dataset and layer have to live in the same context,
if dataset is deleted all layer references will get lost
.. versionadded:: 0.12.0
Parameters
----------
filename : string
vector file name
driver : string
gdal driver string
Returns
-------
dataset : gdal.Dataset
dataset
layer : ogr.Layer
layer
"""
dataset = gdal.OpenEx(filename)
if driver:
gdal.GetDriverByName(driver)
layer = dataset.GetLayer()
return dataset, layer
def open_shape(filename, driver=None):
"""Open shapefile, return gdal.Dataset and OGR.Layer
.. warning:: dataset and layer have to live in the same context,
if dataset is deleted all layer references will get lost
.. versionadded:: 0.6.0
Parameters
----------
filename : string
shapefile name
driver : string
gdal driver string
Returns
-------
dataset : gdal.Dataset
dataset
layer : ogr.Layer
layer
"""
if driver is None:
driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open(filename)
if dataset is None:
print('Could not open file')
raise IOError
layer = dataset.GetLayer()
return dataset, layer
def open_raster(filename, driver=None):
"""Open raster file, return gdal.Dataset
.. versionadded:: 0.6.0
Parameters
----------
filename : string
raster file name
driver : string
gdal driver string
Returns
-------
dataset : gdal.Dataset
dataset
"""
dataset = gdal.OpenEx(filename)
if driver:
gdal.GetDriverByName(driver)
return dataset
def read_safnwc(filename):
"""Read MSG SAFNWC hdf5 file into a gdal georeferenced object
Parameters
----------
filename : string
satellite file name
Returns
-------
ds : gdal.DataSet
with satellite data
"""
root = gdal.Open(filename)
ds1 = gdal.Open('HDF5:' + filename + '://CT')
ds = gdal.GetDriverByName('MEM').CreateCopy('out', ds1, 0)
# name = os.path.basename(filename)[7:11]
try:
proj = osr.SpatialReference()
proj.ImportFromProj4(ds.GetMetadata()["PROJECTION"])
except Exception:
raise NameError("No metadata for satellite file %s" % filename)
geotransform = root.GetMetadata()["GEOTRANSFORM_GDAL_TABLE"].split(",")
geotransform[0] = root.GetMetadata()["XGEO_UP_LEFT"]
geotransform[3] = root.GetMetadata()["YGEO_UP_LEFT"]
ds.SetProjection(proj.ExportToWkt())
ds.SetGeoTransform([float(x) for x in geotransform])
return ds
def reprojectPanBand(src,match_ds,dst_filename):
src_proj = src.GetProjection()
src_geotrans = src.GetGeoTransform()
match_proj = match_ds.GetProjection()
match_geotrans = match_ds.GetGeoTransform()
wide = match_ds.RasterXSize
high = match_ds.RasterYSize
dst = gdal.GetDriverByName('GTiff').Create(dst_filename, wide, high, 1, gdalconst.GDT_Int16)
dst.SetGeoTransform( match_geotrans )
dst.SetProjection( match_proj)
gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Bilinear)
return(gdal.Open(dst_filename,gdalconst.GA_ReadOnly))
def dump_raster(rast, rast_path, xoff=0, yoff=0, driver='GTiff', nodata=None):
'''
Creates a raster file from a given GDAL dataset (raster). Arguments:
rast A gdal.Dataset; does NOT accept NumPy array
rast_path The path of the output raster file
xoff Offset in the x-direction; should be provided when clipped
yoff Offset in the y-direction; should be provided when clipped
driver The name of the GDAL driver to use (determines file type)
nodata The NoData value; defaults to -9999.
'''
driver = gdal.GetDriverByName(driver)
sink = driver.Create(rast_path, rast.RasterXSize, rast.RasterYSize,
rast.RasterCount, rast.GetRasterBand(1).DataType)
assert sink is not None, 'Cannot create dataset; there may be a problem with the output path you specified'
sink.SetGeoTransform(rast.GetGeoTransform())
sink.SetProjection(rast.GetProjection())
for b in range(1, rast.RasterCount + 1):
dat = rast.GetRasterBand(b).ReadAsArray()
sink.GetRasterBand(b).WriteArray(dat)
sink.GetRasterBand(b).SetStatistics(*map(np.float64,
[dat.min(), dat.max(), dat.mean(), dat.std()]))
if nodata is None:
nodata = rast.GetRasterBand(b).GetNoDataValue()
if nodata is None:
nodata = -9999
sink.GetRasterBand(b).SetNoDataValue(np.float64(nodata))
sink.FlushCache()
def convertLabelStringToPoly(shapeFileSrc, outGeoJSon, labelType='Airplane'):
shapeSrc = ogr.Open(shapeFileSrc)
source_layer = shapeSrc.GetLayer()
source_srs = source_layer.GetSpatialRef()
# Create the output Layer
outDriver = ogr.GetDriverByName("geojson")
if os.path.exists(outGeoJSon):
outDriver.DeleteDataSource(outGeoJSon)
outDataSource = outDriver.CreateDataSource(outGeoJSon)
outLayer = outDataSource.CreateLayer("groundTruth", source_srs, geom_type=ogr.wkbPolygon)
# Add input Layer Fields to the output Layer
inLayerDefn = source_layer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
fieldDefn = inLayerDefn.GetFieldDefn(i)
outLayer.CreateField(fieldDefn)
outLayer.CreateField(ogr.FieldDefn("Length_m", ogr.OFTReal))
outLayer.CreateField(ogr.FieldDefn("Width_m", ogr.OFTReal))
outLayer.CreateField(ogr.FieldDefn("Aspect(L/W)", ogr.OFTReal))
outLayer.CreateField(ogr.FieldDefn("compassDeg", ogr.OFTReal))
outLayerDefn = outLayer.GetLayerDefn()
for inFeature in source_layer:
outFeature = ogr.Feature(outLayerDefn)
for i in range(0, inLayerDefn.GetFieldCount()):
outFeature.SetField(inLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
geom = inFeature.GetGeometryRef()
if labelType == 'Airplane':
poly, Length, Width, Aspect, Direction = evaluateLineStringPlane(geom, label='Airplane')
elif labelType == 'Boat':
poly, Length, Width, Aspect, Direction = evaluateLineStringBoat(geom, label='Boat')
outFeature.SetGeometry(poly)
outFeature.SetField("Length_m", Length)
outFeature.SetField("Width_m", Width)
outFeature.SetField("Aspect(L/W)", Aspect)
outFeature.SetField("compassDeg", Direction)
outLayer.CreateFeature(outFeature)
def reproject_tiff_custom(old_name, new_name, new_x_size, new_y_size,
new_geo_transform, new_projection, unit, mode):
"""
Reprojects an tiff with custom tiff arguments. Can be used to warp tiff.
If no projection is provided, fallback to old projection.
Keyword arguments:
old_name -- the name of the old tiff file
new_name -- the name of the output tiff file
new_x_size -- the number of new size in x dimension
new_y_size -- the number of new size in y dimension
new_geo_transform -- the new geo transform to apply
new_projection -- the new projection to use
unit -- the gdal unit in which the operation will be performed
mode -- the gdal mode used for warping
"""
mem_drv = gdal.GetDriverByName("MEM")
old = gdal.Open(old_name)
r = old.GetRasterBand(1)
r.GetNoDataValue()
# Adds 1 to keep the original zeros (reprojectImage maps NO_DATA to 0)
old_array = old.ReadAsArray()
mask = old_array == old.GetRasterBand(1).GetNoDataValue()
old_array += 1
old_array[mask] = 0
temp = mem_drv.Create("temp", old.RasterXSize, old.RasterYSize, 1, unit)
temp.SetGeoTransform(old.GetGeoTransform())
temp.SetProjection(old.GetProjection())
temp.GetRasterBand(1).WriteArray(old_array)
new = mem_drv.Create(new_name, new_x_size, new_y_size, 1, unit)
new.SetGeoTransform(new_geo_transform)
if new_projection is None:
new.SetProjection(old.GetProjection())
else:
new.SetProjection(new_projection)
res = gdal.ReprojectImage(temp, new, old.GetProjection(), new_projection,
mode)
assert res == 0, 'Error in ReprojectImage'
arr = new.ReadAsArray()
mask = arr != 0
arr -= 1
arr[~mask] = 0
new = None
temp = None
return arr, mask
def makeDensityRaster(speciesOcc, inVector, pixelSize, outRas, noData):
srcVector = ogr.Open(inVector)
srcLayer = srcVector.GetLayer()
srs = srcLayer.GetSpatialRef()
# if the layer is not wgs84
if srs.GetAttrValue("AUTHORITY", 1) != '4326':
print('Layer projection should be WGS84!')
return
xMin, xMax, yMin, yMax = srcLayer.GetExtent()
# Create the destination data source
xRes = int((xMax - xMin) / pixelSize)
yRes = int((yMax - yMin) / pixelSize)
targetRas = gdal.GetDriverByName('GTiff').Create(outRas, xRes, yRes, 1, 6) # 6 == float
targetRas.SetGeoTransform((xMin, pixelSize, 0, yMax, 0, -pixelSize))
band = targetRas.GetRasterBand(1)
band.SetNoDataValue(noData)
# Rasterize clips the raster band
gdal.RasterizeLayer(targetRas, [1], srcLayer, None, None, [0], ['ALL_TOUCHED=TRUE'])
#os.remove(outRas)
g = band.ReadAsArray()
for point in speciesOcc.iterrows():
xi = int((point[1]['x'] - xMin) / pixelSize)
yi = int((point[1]['y'] - yMax) / -pixelSize)
try:
if g[yi,xi] != noData:
g[yi,xi] += 1
else:
print('point ({}, {}) out of bounds'.format(point[1]['x'], point[1]['y']))
except:
print('point ({}, {}) out of bounds'.format(point[1]['x'], point[1]['y']))
pass
band.WriteArray(g)
targetRasSRS = osr.SpatialReference()
targetRasSRS.ImportFromWkt(srs.ExportToWkt())
targetRas.SetProjection(targetRasSRS.ExportToWkt())
band.FlushCache()
def createRaster(outRas, extCells, pixelSize, cellValues = 'lat', dataType = "float32", noData = -9999, inVector = None):
NP2GDAL_CONVERSION = { "uint8": 1, "uint16": 2, "int16": 3,
"uint32": 4, "int32": 5, "float32": 6, "float64": 7,
"complex64": 10, "complex128": 11,
}
if os.path.exists(outRas):
print('Raster file already excists!')
return
# Create the destination data source
xRes = int(numpy.ceil((extCells[2] - extCells[0]) / pixelSize))
yRes = int(numpy.ceil((extCells[3] - extCells[1]) / pixelSize))
targetRas = gdal.GetDriverByName('GTiff').Create(outRas, xRes, yRes, 1, NP2GDAL_CONVERSION[dataType])
targetRas.SetGeoTransform((extCells[0], pixelSize, 0, extCells[3], 0, -pixelSize))
band = targetRas.GetRasterBand(1)
band.SetNoDataValue(noData)
if inVector != None:
srcVector = ogr.Open(inVector)
srcLayer = srcVector.GetLayer()
srs = srcLayer.GetSpatialRef()
# if the layer is not wgs84
if srs.GetAttrValue("AUTHORITY", 1) != '4326':
print('Layer projection should be WGS84!')
os.remove(outRas)
return
# Rasterize clips the raster band
gdal.RasterizeLayer(targetRas, [1], srcLayer, None, None, [0], ['ALL_TOUCHED=TRUE'])
g = band.ReadAsArray()
else:
g = numpy.zeros((yRes,xRes), eval('numpy.{}'.format(dataType)))
#populate matrix with random numbers
for i in range(yRes):
for j in range(xRes):
if g[i,j] != noData:
if cellValues == 'lat':
g[i,j] = i
elif cellValues == 'lon':
g[i,j] = j
elif cellValues == 'random':
g[i,j] = numpy.random.randint(50)
band.WriteArray(g)
targetRasSRS = osr.SpatialReference()
targetRasSRS.ImportFromEPSG(4326)
targetRas.SetProjection(targetRasSRS.ExportToWkt())
band.FlushCache()
#function to filter raster cells based on the coverage by some vector features
def geom2shp(geom, out_fn, fields=False):
"""Write out a new shapefile for input geometry
"""
from pygeotools.lib import timelib
driverName = "ESRI Shapefile"
drv = ogr.GetDriverByName(driverName)
if os.path.exists(out_fn):
drv.DeleteDataSource(out_fn)
out_ds = drv.CreateDataSource(out_fn)
out_lyrname = os.path.splitext(os.path.split(out_fn)[1])[0]
geom_srs = geom.GetSpatialReference()
geom_type = geom.GetGeometryType()
out_lyr = out_ds.CreateLayer(out_lyrname, geom_srs, geom_type)
if fields:
field_defn = ogr.FieldDefn("name", ogr.OFTString)
field_defn.SetWidth(128)
out_lyr.CreateField(field_defn)
field_defn = ogr.FieldDefn("path", ogr.OFTString)
field_defn.SetWidth(254)
out_lyr.CreateField(field_defn)
#field_defn = ogr.FieldDefn("date", ogr.OFTString)
#This allows sorting by date
field_defn = ogr.FieldDefn("date", ogr.OFTInteger)
field_defn.SetWidth(32)
out_lyr.CreateField(field_defn)
field_defn = ogr.FieldDefn("decyear", ogr.OFTReal)
field_defn.SetPrecision(8)
field_defn.SetWidth(64)
out_lyr.CreateField(field_defn)
out_feat = ogr.Feature(out_lyr.GetLayerDefn())
out_feat.SetGeometry(geom)
if fields:
#Hack to force output extesion to tif, since out_fn is shp
out_path = os.path.splitext(out_fn)[0] + '.tif'
out_feat.SetField("name", os.path.split(out_path)[-1])
out_feat.SetField("path", out_path)
#Try to extract a date from input raster fn
out_feat_date = timelib.fn_getdatetime(out_fn)
if out_feat_date is not None:
datestamp = int(out_feat_date.strftime('%Y%m%d'))
#out_feat_date = int(out_feat_date.strftime('%Y%m%d%H%M'))
out_feat.SetField("date", datestamp)
decyear = timelib.dt2decyear(out_feat_date)
out_feat.SetField("decyear", decyear)
out_lyr.CreateFeature(out_feat)
out_ds = None
#return status?
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 creaLayer(self):
layer = self.iface.activeLayer()
provider = layer.dataProvider()
path = provider.dataSourceUri()
(raiz, filename) = os.path.split(path)
dataset = gdal.Open(path)
# Get projection
prj = dataset.GetProjection()
# setting band
number_band = 1
band = dataset.GetRasterBand(number_band)
# Get raster metadata
geotransform = dataset.GetGeoTransform()
# Set name of output raster
output_file = "C:\\Users\\caligola\\Desktop\\s\\raster_output.tif"
# Create gtif file with rows and columns from parent raster
driver = gdal.GetDriverByName("GTiff")
raster = self.changeRasterValues(band)
dst_ds = driver.Create(output_file,
band.XSize,
band.YSize,
number_band,
band.DataType)
# writting output raster
dst_ds.GetRasterBand(number_band).WriteArray(raster)
# setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(geotransform)
# setting spatial reference of output raster
srs = osr.SpatialReference(wkt=prj)
dst_ds.SetProjection(srs.ExportToWkt())
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)
