processor.py 6.15 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#!/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
25
26
from rasterio import features
import geopandas
27
28
29
30
31
32
33
34
35
36

import os

# Initialize log
logger = logging.getLogger(__name__)


class TileProcessor:
    @staticmethod
    def intersect_datasource_with_tile(datasource, tile):
37
        """Returns a subset of the datasource.raster_files_index with the raster tiles
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
        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)

58
            >>> intersect_datasource_with_tile(datasource, tile)
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

        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):
75
        """Returns the subset of the raster and its affine transformation based
76
77
78
79
80
81
82
83
        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:
84
            pixel_values (list of numpy.ndarray): List of arrays with pixel values that
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
                                    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
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161

    @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