Commit 72580b3a authored by pillot's avatar pillot
Browse files

Initial Commit

parents
The MIT Licence
Copyright (c) 2020 Benjamin Pillot
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NON INFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
# -*- coding: utf-8 -*-
""" Module summary description.
More detailed description.
"""
__version__ = '1.0'
__author__ = 'Benjamin Pillot'
__copyright__ = 'Copyright 2020, Benjamin Pillot'
__email__ = 'benjaminpillot@riseup.net'
import gdal
import osr
gdal.UseExceptions()
osr.UseExceptions()
# -*- coding: utf-8 -*-
""" Module summary description.
More detailed description.
"""
from functools import wraps
from pyraster.tools.conversion import _resample_raster, _padding
from pyraster.tools.merge import _merge
from pyraster.utils import lazyproperty
def return_new_instance(method):
@wraps(method)
def _return_new_instance(self, *args, **kwargs):
output = method(self, *args, **kwargs)
try:
new_self = self.__class__(output)
except RuntimeError:
return output
return _return_new_instance
class RasterBase:
crs = None
file = None
gdal_dataset = None
gdal_driver = None
rasterio_dataset = None
def __init__(self):
pass
@classmethod
def merge(cls, rasters, bounds=None):
return _merge(rasters, bounds)
@return_new_instance
def pad_extent(self, pad_x, pad_y, value):
""" Pad raster extent with given values
Description
-----------
Pad raster extent, i.e. add pad value around raster bounds
Parameters
----------
pad_x: int
x padding size (new width will therefore be RasterXSize + 2 * pad_x)
pad_y: int
y padding size (new height will therefore be RasterYSize + 2 * pad_y)
value: int or float
value to set to pad area around raster
Returns
-------
RasterBase:
New instance
"""
return _padding(self, pad_x, pad_y, value)
@return_new_instance
def resample(self, factor):
""" Resample raster
Description
-----------
Resample raster with respect to resampling factor.
The higher the factor, the higher the resampling.
Parameters
----------
factor: int or float
Resampling factor
Returns
-------
RasterBase:
New temporary resampled instance
"""
return _resample_raster(self, factor)
@lazyproperty
def bounds(self):
return self.rasterio_dataset.bounds
@lazyproperty
def nb_band(self):
return self.gdal_dataset.RasterCount
@lazyproperty
def no_data(self):
return [self.gdal_dataset.GetRasterBand(band + 1).GetNoDataValue() for band in range(self.nb_band)]
@lazyproperty
def resolution(self):
return self.rasterio_dataset.res
@lazyproperty
def x_origin(self):
return self.rasterio_dataset.bounds[0]
@lazyproperty
def x_size(self):
return self.gdal_dataset.RasterXSize
@lazyproperty
def y_origin(self):
return self.rasterio_dataset.bounds[3]
@lazyproperty
def y_size(self):
return self.gdal_dataset.RasterYSize
# -*- coding: utf-8 -*-
""" Functions related to CRS conversion and computation
More detailed description.
"""
import osr
import pyproj
def is_equal_proj(proj1, proj2):
""" Compare 2 projections
Parameters
----------
proj1: int or str or dict or pyproj.Proj
valid projection name
proj2: int or str or dict or pyproj.Proj
valid projection name
Returns
-------
boolean:
True or False
"""
# From an idea from https://github.com/jswhit/pyproj/issues/15
# Use OGR library to compare projections
srs = [srs_from(proj1), srs_from(proj2)]
return bool(srs[0].IsSame(srs[1]))
def proj4_from_wkt(wkt):
""" Convert wkt srs to proj4
Description
-----------
Parameters
----------
wkt:
Returns
-------
"""
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
return srs.ExportToProj4()
def proj4_from(proj):
""" Convert projection to proj4 string
Description
-----------
Convert projection string, dictionary, etc.
to proj4 string
Parameters
----------
proj:
Returns
-------
"""
if type(proj) == int:
try:
proj4_str = pyproj.Proj('epsg:%d' % proj).srs
except (ValueError, RuntimeError):
raise ValueError("Invalid EPSG code")
elif type(proj) == str or type(proj) == dict:
try:
proj4_str = pyproj.Proj(proj).srs
except RuntimeError:
try:
proj4_str = proj4_from_wkt(proj)
except (RuntimeError, TypeError):
raise ValueError("Invalid projection string or dictionary")
elif type(proj) == pyproj.Proj:
proj4_str = proj.srs
else:
raise ValueError("Invalid projection format: '{}'".format(type(proj)))
return proj4_str
def srs_from(proj):
""" Get spatial reference system from projection
Description
-----------
Parameters
----------
proj:
Returns
-------
SpatialReference instance (osgeo.osr package)
"""
proj4 = proj4_from(proj)
srs = osr.SpatialReference()
srs.ImportFromProj4(proj4)
return srs
def wkt_from(proj):
""" Get WKT spatial reference system from projection
Description
-----------
Parameters
----------
proj:
Returns
-------
"""
return srs_from(proj).ExportToWkt()
# -*- coding: utf-8 -*-
""" Module summary description.
More detailed description.
"""
import os
from tempfile import mkstemp
class File:
def __init__(self, path):
self.path = path
def __enter__(self):
return self.path
class TempFile(File):
def __del__(self):
try:
os.remove(self.path)
except FileNotFoundError:
pass
class MemmapFile(TempFile):
def __init__(self):
super().__init__(mkstemp(suffix='.dat')[1])
def __exit__(self, exc_type, exc_val, exc_tb):
del self
class RasterTempFile(TempFile):
""" Create temporary raster file
"""
def __init__(self):
super().__init__(mkstemp(suffix='.tif')[1])
# -*- coding: utf-8 -*-
""" Module summary description.
More detailed description.
"""
# -*- coding: utf-8 -*-
""" Module summary description.
More detailed description.
"""
import rasterio
from pyraster.tools.merge import _rasterio_merge_modified
r1 = rasterio.open("/home/benjamin/Documents/pyraster/pyraster/test/raster_1.tif")
r2 = rasterio.open("/home/benjamin/Documents/pyraster/pyraster/test/raster_2.tif")
sources = [r1, r2]
out_file = "/home/benjamin/Documents/pyraster/pyraster/test/test_merge.tif"
_rasterio_merge_modified(sources, out_file)
# -*- coding: utf-8 -*-
""" Module summary description.
More detailed description.
"""
__author__ = 'Benjamin Pillot'
__copyright__ = 'Copyright 2020, Benjamin Pillot'
__email__ = 'benjaminpillot@riseup.net'
def set_no_data(gdal_ds, no_data):
""" Set no data value into gdal dataset
Description
-----------
Parameters
----------
gdal_ds: gdal.Dataset
gdal dataset
no_data: list or tuple
list of no data values corresponding to each raster band
"""
for band in range(gdal_ds.RasterCount):
try:
gdal_ds.GetRasterBand(band + 1).SetNoDataValue(no_data[band])
except TypeError:
pass
# -*- coding: utf-8 -*-
""" Module summary description.
More detailed description.
"""
import gdal
from pyraster.crs import srs_from
from pyraster.io import RasterTempFile
from pyraster.tools import set_no_data
def _padding(raster, pad_x, pad_y, pad_value):
""" Add pad values around raster
Description
-----------
Parameters
----------
raster: RasterBase
raster to pad
pad_x: int
x padding size (new width will therefore be RasterXSize + 2 * pad_x)
pad_y: int
y padding size (new height will therefore be RasterYSize + 2 * pad_y)
pad_value: int or float
value to set to pad area around raster
Returns
-------
"""
with RasterTempFile() as out_file:
out_ds = raster.gdal_driver.Create(out_file, raster.x_size + 2 * pad_x, raster.y_size + 2 * pad_y,
raster.nb_band, raster.gdal_dataset.GetRasterBand(1).DataType)
out_ds.SetGeoTransform((raster.x_origin - pad_x * raster.resolution[0], raster.resolution[0],
0, raster.y_origin + pad_y * raster.resolution[1], 0, -raster.resolution[1]))
out_ds.SetProjection(raster.gdal_dataset.GetProjection())
set_no_data(out_ds, raster.no_data)
for band in range(1, raster.nb_band + 1):
out_ds.GetRasterBand(band).Fill(pad_value)
gdal.Warp(out_ds, raster.gdal_dataset)
# Close dataset
out_ds = None
return out_file
def _project_raster(raster, new_crs):
""" Project raster onto new CRS
"""
with RasterTempFile() as out_file:
gdal.Warp(out_file, raster.gdal_dataset, dstSRS=srs_from(new_crs))
def _resample_raster(raster, factor):
""" Resample raster
Parameters
----------
raster: RasterBase
raster to resample
factor: int or float
Resampling factor
"""
with RasterTempFile() as out_file:
out_ds = raster.gdal_driver.Create(out_file, raster.x_size * factor, raster.y_size * factor, raster.nb_band,
raster.gdal_dataset.GetRasterBand(1).DataType)
out_ds.SetGeoTransform((raster.x_origin, raster.resolution[0] / factor, 0, raster.y_origin, 0,
-raster.resolution[1] / factor))
out_ds.SetProjection(raster.gdal_dataset.GetProjection())
set_no_data(out_ds, raster.no_data)
for band in range(1, raster.nb_band+1):
gdal.RegenerateOverview(raster.gdal_dataset.GetRasterBand(band), out_ds.GetRasterBand(band), 'mode')
# Close dataset
out_ds = None
return out_file
# -*- coding: utf-8 -*-
""" Module summary description.
More detailed description.
"""
import numpy as np
import rasterio
from rasterio import windows
from rasterio.transform import Affine
def _merge(rasters, bounds):
pass
def _rasterio_merge_modified(sources, out_file, bounds=None, driver="GTiff", precision=7):
""" Modified rasterio merge
Description
-----------
Merge set of rasters using modified
rasterio merging tool. Only import as
numpy arrays the source datasets, and
write each source to destination raster.
Parameters
----------
sources: list
list of rasterio datasets
out_file: str
valid path to the raster file to be written
bounds: tuple
valid boundary tuple (optional)
driver: str
valid gdal driver (optional)
precision: int
float precision (optional)
Note
----
Adapted from
https://gis.stackexchange.com/questions/348925/merging-rasters-with-rasterio-in-blocks-to-avoid-memoryerror
"""
# adapted from https://github.com/mapbox/rasterio/blob/master/rasterio/merge.py
first = sources[0]
first_res = first.res
dtype = first.dtypes[0]
# Determine output band count
output_count = first.count
# Extent of all inputs
if bounds:
dst_w, dst_s, dst_e, dst_n = bounds
else:
# scan input files
xs = []
ys = []
for src in sources:
left, bottom, right, top = src.bounds
xs.extend([left, right])
ys.extend([bottom, top])
dst_w, dst_s, dst_e, dst_n = min(xs), min(ys), max(xs), max(ys)
out_transform = Affine.translation(dst_w, dst_n)
# Resolution/pixel size
res = first_res
out_transform *= Affine.scale(res[0], -res[1])
# Compute output array shape. We guarantee it will cover the output
# bounds completely
output_width = int(np.ceil((dst_e - dst_w) / res[0]))
output_height = int(np.ceil((dst_n - dst_s) / res[1]))
# Adjust bounds to fit
dst_e, dst_s = out_transform * (output_width, output_height)
# create destination array
# destination array shape
shape = (output_height, output_width)
dest_profile = {
"driver": driver,
"height": shape[0],
"width": shape[1],
"count": output_count,
"dtype": dtype,
"crs": sources[0].crs.to_proj4(),
"nodata": sources[0].nodata,
"transform": out_transform
}
# open output file in write/read mode and fill with destination mosaick array
with rasterio.open(out_file, 'w+', **dest_profile) as mosaic_raster:
for src in sources:
# 1. Compute spatial intersection of destination and source
src_w, src_s, src_e, src_n = src.bounds
int_w = src_w if src_w > dst_w else dst_w
int_s = src_s if src_s > dst_s else dst_s
int_e = src_e if src_e < dst_e else dst_e
int_n = src_n if src_n < dst_n else dst_n
# 2. Compute the source window
src_window = windows.from_bounds(
int_w, int_s, int_e, int_n, src.transform, precision=precision)
src_window = src_window.round_shape()
# 3. Compute the destination window
dst_window = windows.from_bounds(
int_w, int_s, int_e, int_n, out_transform, precision=precision)
dst_window = windows.Window(int(round(dst_window.col_off)),
int(round(dst_window.row_off)),
int(round(dst_window.width)),
int(round(dst_window.height)))
out_shape = (dst_window.height, dst_window.width)
for band in range(1, output_count+1):
src_array = src.read(band, out_shape=out_shape, window=src_window)
# before writing the window, replace source nodata with dest nodata as it can
# already have been written (e.g. another adjacent country)
dst_array = mosaic_raster.read(band, out_shape=out_shape, window=dst_window)
mask = src_array == src.nodata
src_array[mask] = dst_array[mask]
mosaic_raster.write(src_array, band, window=dst_window)
# -*- coding: utf-8 -*-
""" Module summary description.
More detailed description.
"""
def lazyproperty(func):
name = '_lazy_' + func.__name__
@property
def lazy(self):
if hasattr(self, name):
return getattr(self, name)
else:
value = func(self)
setattr(self, name, value)
return value
return lazy