import glob
import logging
import os
import numpy as np
from osgeo import gdal, osr
from ..var_config import nan_value
[docs]
def open_raster(file_name, band_number=1):
"""Opens a raster file and accesses its bands.
Args:
file_name (str): The raster file directory and name.
band_number (int): The Raster band number to open (default: ``1``).
Returns:
osgeo.gdal.Dataset: A raster dataset a Python object.
osgeo.gdal.Band: The defined raster band as Python object.
"""
gdal.UseExceptions()
# open raster file or return None if not accessible
try:
raster = gdal.Open(file_name)
except RuntimeError as e:
logging.error("Cannot open raster.")
print(e)
return nan_value, nan_value
# open raster band or return None if corrupted
try:
raster_band = raster.GetRasterBand(band_number)
except RuntimeError as e:
logging.error("Cannot access raster band.")
logging.error(e)
return raster, nan_value
return raster, raster_band
[docs]
def create_raster(
file_name,
raster_array,
bands=1,
origin=None,
epsg=4326,
pixel_width=10.0,
pixel_height=10.0,
nan_val=nan_value,
rdtype=gdal.GDT_Float32,
geo_info=False,
rotation_angle=None,
shear_pixels=True,
options=None,
):
"""Converts an ``ndarray`` (``numpy.array``) to a GeoTIFF raster.
Args:
file_name (str): Target file name, including directory; must end on ``".tif"``.
raster_array (``ndarray`` or ``list``): 2d-array (no bands) or list (bands) of 2-darrays of values to rasterize. If a list of 2d-arrays is provided, the length of the list will correspond to the number of bands added to the raster (supersedes ``bands``).
bands (int): Number of bands to write to the raster (default: ``1``).
origin (tuple): Coordinates (x, y) of the origin.
epsg (int): EPSG:XXXX projection to use (default: ``4326``).
pixel_height (float OR int): Pixel height as multiple of the base units defined with the EPSG number (default: ``10`` meters).
pixel_width (float OR int): Pixel width as multiple of the base units defined with the EPSG number (default: ``10`` meters).
nan_val (``int`` or ``float``): No-data value to be used in the raster. Replaces non-numeric and ``np.nan`` in the ``ndarray``. (default: ``geoconfig.nan_value``).
rdtype: `gdal.GDALDataType <https://gdal.org/doxygen/gdal_8h.html#a22e22ce0a55036a96f652765793fb7a4>`_ raster data type (default: gdal.GDT_Float32 (32 bit floating point).
geo_info (tuple): Defines a ``gdal.DataSet.GetGeoTransform`` object and supersedes ``origin``, ``pixel_width``, ``pixel_height`` (default: ``False``).
rotation_angle (float): Rotate (in degrees) not North-up rasters. The default value (``0``) corresponds to north-up (only modify if you know what you are doing).
shear_pixels (bool): Use with ``rotation_angle`` to shear pixels as well (default: ``True``).
options (list): Raster creation options - default is ['PROFILE=GeoTIFF']. Add 'PHOTOMETRIC=RGB' to create an RGB image raster.
Returns:
int: ``0`` if successful, otherwise ``-1``.
Hint:
For processing airborne imagery, the ``rotation_angle`` corresponds to the bearing angle of the aircraft with reference to true, not magnetic North.
"""
gdal.UseExceptions()
if options is None:
options = ["PROFILE=GeoTIFF"]
# check out driver
driver = gdal.GetDriverByName("GTiff")
# create raster dataset with number of cols and rows of the input array
try:
# overwrite number of bands if multiple arrays are provided in a list
if isinstance(raster_array, list):
bands = len(raster_array)
cols = raster_array[0].shape[1]
rows = raster_array[0].shape[0]
else:
cols = raster_array.shape[1]
rows = raster_array.shape[0]
except TypeError:
logging.error("Provided array is not a numpy.ndarray.")
return -1
try:
logging.info(" * creating new raster with %1i bands ...", bands)
new_raster = driver.Create(file_name, cols, rows, bands, eType=rdtype, options=options)
except RuntimeError:
logging.error("Could not create %s.", str(file_name))
return -1
# apply geo-origin and pixel dimensions
if not geo_info:
try:
origin_x = origin[0]
origin_y = origin[1]
except IndexError:
logging.error(
"Wrong origin format (required: (INT, INT) - provided: %s).", str(origin)
)
return -1
if rotation_angle:
try:
logging.info(" * rotating image by %0.2f deg", float(rotation_angle))
except ValueError:
logging.error(
"The provided rotation angle is not a number. Re-try with a numeric rotation angle (in degrees)."
)
return -1
rotation_angle = np.deg2rad(rotation_angle)
x_rotation = -1 * pixel_width * np.sin(rotation_angle)
y_rotation = -pixel_height * np.cos(rotation_angle)
if shear_pixels:
pixel_width = pixel_width * np.cos(rotation_angle)
pixel_height = pixel_height * np.sin(rotation_angle)
else:
x_rotation = 0.0
y_rotation = 0.0
try:
new_raster.SetGeoTransform(
(origin_x, pixel_width, x_rotation, origin_y, y_rotation, -pixel_height)
)
except RuntimeError:
logging.error(
"Invalid origin (must be INT) or pixel_height/pixel_width (must be INT) provided."
)
return -1
else:
try:
new_raster.SetGeoTransform(geo_info)
except RuntimeError as e:
logging.error(e)
return -1
# write array contents to band(s)
for b in range(bands):
if isinstance(raster_array, list):
# use array item of list if multiple arrays provided
write_array = raster_array[b]
else:
write_array = raster_array
# replace np.nan values
write_array[np.isnan(write_array)] = nan_val
band = new_raster.GetRasterBand(b + 1)
band.SetNoDataValue(nan_val)
band.WriteArray(write_array)
band.SetScale(1.0)
# release band
band.FlushCache()
# create projection and assign to raster
srs = osr.SpatialReference()
try:
srs.ImportFromEPSG(epsg)
except RuntimeError as e:
logging.error(e)
return -1
new_raster.SetProjection(srs.ExportToWkt())
logging.info(" * successfully created %s", file_name)
return 0
[docs]
def xy_raster_shift(
file_name,
x_shift,
y_shift,
bands=1,
rdtype=gdal.GDT_Float32,
nan_val=nan_value,
compress=True,
options=None,
compress_config=None,
):
"""Creates new geotiff raster with shifts in x and y direction. If enabled compresses it also compresses file.
Args:
file_name (string): File path of GeoTiff
x_shift (float OR int): Shift origin in x direction. *Check that correct units are used. Example: wgs 84 is in degrees
y_shift (float OR int): Shift origin in y direction. *Check that correct units are used. Example: wgs 84 is in degrees
bands (int): Number of bands default is 1, however check raster to see how many are required. Example: RGBA=4
rdtype: `gdal.GDALDataType <https://gdal.org/doxygen/gdal_8h.html#a22e22ce0a55036a96f652765793fb7a4>`_ raster data type (default: gdal.GDT_Float32 (32 bit floating point).
nan_val (``int`` or ``float``): No-data value to be used in the raster. Replaces non-numeric and ``np.nan`` in the ``ndarray``. (default: ``geoconfig.nan_value``).
compress (Bool): If True creates compressed version of the GeoTiff
options (list): Raster creation options - default is ['PROFILE=GeoTIFF']. Add 'PHOTOMETRIC=RGB' to create an RGB image raster.
compress_config: (list) Compress creation options - default is ["COMPRESS=LZW", "TILED=YES"] LZW=Lempel-Ziv-Welch-Algorithm See gdal.Translate for more options
Returns:
int: ``0`` if successful, otherwise ``-1``.
Hint:
For drone rasters try
bands=4 (rgba) rdtype=gdal.GDT_Byte nan_val=0 options=['PROFILE=GeoTIFF','PHOTOMETRIC=RGB'])
Bugs: Issues displaying logging
"""
if options is None:
options = ["PROFILE=GeoTIFF"]
if compress_config is None:
compress_config = ["COMPRESS=LZW", "TILED=YES"]
# Gdal opens Tiff
try:
tif = gdal.Open(file_name)
except RuntimeError:
logging.error("Cannot open raster in Gdal.")
return -1
# Extracting geo_transform
geo_transform = tif.GetGeoTransform()
# Extracting arrays from raster to create an array list
try:
list_array = []
for n in range(bands):
list_array.append(tif.GetRasterBand(n + 1).ReadAsArray())
except RuntimeError:
logging.error("Cannot create array list from bands.")
return -1
# sets the geodata into individual types
try:
origin_x = geo_transform[0]
origin_y = geo_transform[3]
pixel_width = geo_transform[1]
pixel_height = geo_transform[5]
pixel_height = (
pixel_height * -1
) # setting negative so raster is not mirrored over y-intercept, due to create_raster's "*-1 pixelheight"
proj = osr.SpatialReference(tif.GetProjection())
epsg = int(proj.GetAttrValue("AUTHORITY", 1))
except ValueError:
logging.error("Problems with geodata")
try:
logging.info(
f"Creating raster with x shift of {float(x_shift)} units and y shift of {float(y_shift)} units"
)
except ValueError:
logging.error("The provided x and y shifts are not numbers.")
return -1
# Gdal Shift
origin = (origin_x + x_shift, origin_y + y_shift)
file_name_new = file_name.replace(".tif", f"shifted x_{x_shift} y_{y_shift}.tif")
try:
create_raster(
file_name_new,
list_array,
bands,
origin,
epsg,
pixel_width,
pixel_height,
nan_val=nan_val,
rdtype=rdtype,
options=options,
)
logging.info("Successfully created shifted raster in " + file_name_new)
except RuntimeError:
logging.error("Could not create raster using geotools.")
return -1
if compress:
try:
logging.info("Creating compressed raster to reduce size")
outfn = file_name_new.replace(".tif", "_compressed.tif")
ds = gdal.Translate(outfn, file_name_new, creationOptions=compress_config)
ds = None
logging.info("Successfully created compressed tiff in " + outfn)
except RuntimeError:
logging.error("Unable to preform compression")
return 0
[docs]
def raster2array(file_name, band_number=1):
"""Extracts a numpy ``ndarray`` from a raster.
Args:
file_name (str): Target file name, including directory; must end on ``".tif"``.
band_number (int): The raster band number to open (default: ``1``).
Returns:
list: three-elements of [``osgeo.DataSet`` of the raster,
``numpy.ndarray`` of the raster ``band_number`` (input) where no-data
values are replaced with ``np.nan``, ``osgeo.GeoTransform`` of
the original raster]
"""
# open the raster and band (see above)
raster, band = open_raster(file_name, band_number=band_number)
try:
# read array data from band
band_array = band.ReadAsArray()
except AttributeError:
logging.error("Could not read array of raster band type=%s.", str(type(band)))
return raster, band, nan_value
try:
# overwrite NoDataValues with np.nan
band_array = np.where(band_array == band.GetNoDataValue(), np.nan, band_array)
except AttributeError:
logging.error("Could not get NoDataValue of raster band type=%s.", str(type(band)))
return raster, band, nan_value
# return the array and GeoTransformation used in the original raster
return raster, band_array, raster.GetGeoTransform()
[docs]
def remove_tif(file_name):
"""Removes a GeoTIFF and its dependent files (e.g., xml).
Args:
file_name (str): Directory (path) and name of a GeoTIFF
Returns:
None: Removes the provided ``file_name`` and all dependencies.
"""
for file in glob.glob(f"{file_name.split('.tif')[0]}*"):
try:
os.remove(file)
except PermissionError:
print(f"WARNING: Could not remove {file} (locked by other program).")
except FileNotFoundError:
print(f"WARNING: The file {file} does not exist.")
[docs]
def clip_raster(polygon, in_raster, out_raster):
"""Clips a raster to a polygon.
Args:
polygon (str): A polygon shapefile name, including directory; must end on ``".shp"``.
in_raster (str): Name of the raster to be clipped, including its directory.
out_raster (str): Name of the target raster, including its directory.
Returns:
None: Creates a new, clipped raster defined with ``out_raster``.
"""
gdal.Warp(out_raster, in_raster, cutlineDSName=polygon)
