#!/usr/bin/env python3
# Copyright (c) 2021:
# Helmholtz-Zentrum Potsdam Deutsches GeoForschungsZentrum GFZ
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero
# General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see http://www.gnu.org/licenses/.
import logging
import rasterio
from rasterio.mask import mask
from shapely.geometry import mapping
from rasterio import features
import geopandas
import os
# Initialize log
logger = logging.getLogger(__name__)
class TileProcessor:
@staticmethod
def intersect_datasource_with_tile(datasource, tile):
"""Returns a subset of the datasource.raster_files_index with the raster tiles
that intersect the given tile.geometry
Args:
datasource (datasource.DataSource): DataSource instance with
crs, pathname and raster_files_index attributes.
tile (tile.Tile): Tile object with quadkey, crs and geometry attributes
Examples:
>>> from obmgapanalysis.datasource import DataSource
>>> from obmgapanalysis.tile import Tile
>>> datasource = DataSource(
crs="epsg:3857",
pathname=os.environ["OBMGAPANALYSIS_BASE_PATH"],
raster_index="GHS_TEST_INDEX.shp"
)
>>> tile = Tile("122100203023320202", datasource.crs)
>>> intersect_datasource_with_tile(datasource, tile)
Returns:
index_list (list): List of the filepaths of the raster files that
intersect with the Tile.geometry.
"""
raster_files_index = datasource.raster_files_index
raster_files_index = raster_files_index[raster_files_index.intersects(tile.geometry)]
raster_files_index["location"] = raster_files_index["location"].apply(
lambda x: os.path.join(datasource.pathname, x)
)
index_list = raster_files_index.location.tolist()
return index_list
@staticmethod
def get_raster_pixels_in_tile(raster_filepaths, tile):
"""Returns the subset of the raster and its affine transformation based
on which pixels intersect with the cell's geometry
Args:
raster_filepaths (list): List of filepaths (may be only one) to raster files
tile (tile.Tile): Tile object with quadkey, crs and geometry attributes
Returns:
pixel_values (list of numpy.ndarray): List of arrays with pixel values that
intersect with the tile.geometry
pixel_georeferences (list of affine.Affine): List with affine transformation
matrices of pixel_values (allows geocoding)
"""
pixel_values = []
pixel_georeferences = []
for i, filepath in enumerate(raster_filepaths):
with rasterio.open(filepath) as src:
feature = [mapping(tile.geometry)]
masked_band = mask(src, feature, crop=True, all_touched=True)
pixel_values.append(masked_band[0].reshape(-1, masked_band[0].shape[-1]))
pixel_georeferences.append(masked_band[1])
return pixel_values, pixel_georeferences
@staticmethod
def polygonize_array(pixel_values, pixel_georeferences, datasource):
"""Returns a Shapely geometry object with built area (multi)polygon for a tile
based on the key pixel_values for the given datasource
Args:
pixel_values (list of numpy.ndarray): List of arrays with pixel values that
intersect with the tile.geometry
pixel_georeferences (list of affine.Affine): List with affine transformation
matrices of pixel_values (allows geocoding)
datasource (datasource.DataSource): DataSource instance with
crs, pathname and raster_files_index attributes.
Returns:
built_area_geometry (shapely.geometry.multipolygon.MultiPolygon): Shapely
polygon or multipolygon
"""
built_pixel_values = datasource.built_pixel_values
results = []
# Run through the list of matrices and extract built up features
for i in range(len(pixel_values)):
features_in_matrix = features.shapes(
pixel_values[i], transform=pixel_georeferences[i]
)
for geometry, pixel_value in features_in_matrix:
if pixel_value in built_pixel_values:
result = {
"properties": {"pixel_value": pixel_value},
"geometry": geometry,
}
results.append(result)
# Completely dissolve all features into one geometry object
if results:
results_geodataframe = geopandas.GeoDataFrame.from_features(results)
built_area_geometry = results_geodataframe.unary_union
return built_area_geometry
@staticmethod
def clip_to_tile_extent(polygon, tile):
"""Returns a Shapely geometry object of a (multi)polygon based on the
intersection of the given geometry (built areas) with the tile.geometry
Args:
polygon (shapely.geometry.multipolygon.MultiPolygon): Shapely polygon
or multipolygon
tile (tile.Tile): Tile object with quadkey, crs and geometry attributes
Returns:
polygon (shapely.geometry.multipolygon.MultiPolygon): Shapely polygon
or multipolygon clipped to the shape of tile.geometry
"""
clipped_polygon = polygon.intersection(tile.geometry)
return clipped_polygon