Commit 88957904 authored by Nicolas Garcia Ospina's avatar Nicolas Garcia Ospina
Browse files

Included processor.py with two basic functions.

parent ea631acd
Pipeline #19348 passed with stage
in 3 minutes and 50 seconds
...@@ -11,7 +11,7 @@ cache: ...@@ -11,7 +11,7 @@ cache:
before_script: before_script:
- apt-get update -y - apt-get update -y
- apt-get install -y python3-pip proj-bin libgeos-dev - apt-get install -y python3-pip proj-bin libgeos-dev gdal-bin libgdal-dev
- python3 -V - python3 -V
- pip3 install virtualenv - pip3 install virtualenv
- rm -rf venv - rm -rf venv
......
#!/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
import os
# Initialize log
logger = logging.getLogger(__name__)
class TileProcessor:
@staticmethod
def intersect_datasource_with_tile(datasource, tile):
"""Return 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)
>>> intersecting_raster_tiles(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):
"""Return 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 matrices 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
...@@ -36,6 +36,7 @@ setup( ...@@ -36,6 +36,7 @@ setup(
"fiona", "fiona",
"rtree", "rtree",
"geopandas", "geopandas",
"rasterio",
], ],
extras_require={ extras_require={
"tests": tests_require, "tests": tests_require,
......
...@@ -22,7 +22,7 @@ import geopandas ...@@ -22,7 +22,7 @@ import geopandas
import os import os
def test_DataSource_init(): def test_init():
datasource = DataSource( datasource = DataSource(
crs="epsg:3857", crs="epsg:3857",
...@@ -31,4 +31,4 @@ def test_DataSource_init(): ...@@ -31,4 +31,4 @@ def test_DataSource_init():
) )
assert datasource.crs == "epsg:3857" assert datasource.crs == "epsg:3857"
isinstance(datasource.raster_files_index, geopandas.GeoDataFrame) assert isinstance(datasource.raster_files_index, geopandas.GeoDataFrame)
#!/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/.
from obmgapanalysis.processor import TileProcessor
from obmgapanalysis.tile import Tile
from obmgapanalysis.datasource import DataSource
import os
import numpy
import affine
# Create an instance for a datasource and a tile
datasource = DataSource(
crs="epsg:3857",
pathname=os.environ["OBMGAPANALYSIS_BASE_PATH"],
raster_files_index="GHS_TEST_INDEX.shp",
)
tile = Tile("122100200320321022", "epsg:3857")
def test_intersect_datasource_with_tile():
# Assert if both classes have the same crs
assert datasource.crs == tile.crs
query = TileProcessor.intersect_datasource_with_tile(datasource=datasource, tile=tile)
# Assert if the resulting raster_files_index subset has 2 entries
assert len(query) == 2
def test_get_raster_pixels_in_tile():
# Find raster files to work with
sample_1 = numpy.array(
[[2, 2, 2, 5], [5, 5, 5, 5], [5, 5, 5, 2], [5, 5, 5, 5], [2, 5, 5, 5], [2, 2, 5, 2]]
)
query = TileProcessor.intersect_datasource_with_tile(datasource=datasource, tile=tile)
# Extract pixels and geocoding information
pixel_values, pixel_georeferences = TileProcessor.get_raster_pixels_in_tile(query, tile)
# Assert datatypes
assert isinstance(pixel_values[0], numpy.ndarray)
assert isinstance(pixel_georeferences[1], affine.Affine)
# Assert content of matrices
assert pixel_values[1].all() == sample_1.all()
assert pixel_georeferences[0][2] == 2549910.0
...@@ -20,7 +20,7 @@ ...@@ -20,7 +20,7 @@
from obmgapanalysis.tile import Tile from obmgapanalysis.tile import Tile
def test_Tile_init(): def test_init():
tile = Tile("122100203310222101", "epsg:3857") tile = Tile("122100203310222101", "epsg:3857")
polygon_coords = ( polygon_coords = (
"POLYGON ((" "POLYGON (("
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment