test_processor.py 4.36 KB
Newer Older
1
2
#!/usr/bin/env python3

3
# Copyright (C) 2021:
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#   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

27

28
29
30
31
32
# 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",
33
    built_pixel_values=[6, 5, 4, 3],
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
)
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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114


def test_poligonize_array():

    query = TileProcessor.intersect_datasource_with_tile(datasource=datasource, tile=tile)
    pixel_values, pixel_georeferences = TileProcessor.get_raster_pixels_in_tile(query, tile)

    expected_geometry = (
        "MULTIPOLYGON (((2549910 4629720, 2549940 4629720, "
        "2549970 4629720, 2549970 4629660, 2549940 4629660, "
        "2549910 4629660, 2549910 4629720)), ((2549970 4629780, "
        "2550000 4629780, 2550090 4629780, 2550090 4629810, "
        "2550120 4629810, 2550120 4629750, 2550090 4629750, "
        "2550090 4629720, 2550120 4629720, 2550120 4629660, "
        "2550090 4629660, 2550090 4629630, 2550060 4629630, "
        "2550060 4629660, 2550030 4629660, 2550030 4629690, "
        "2550000 4629690, 2550000 4629720, 2549970 4629720, "
        "2549970 4629780)))"
    )
    built_geometry = TileProcessor.polygonize_array(
        pixel_values, pixel_georeferences, datasource
    )

    assert str(built_geometry) == expected_geometry


def test_clip_to_tile_extent():

    query = TileProcessor.intersect_datasource_with_tile(datasource=datasource, tile=tile)
    pixel_values, pixel_georeferences = TileProcessor.get_raster_pixels_in_tile(query, tile)

    expected_clipped_geometry = (
        "MULTIPOLYGON (((2549939.26359348 4629720, 2549940 4629720, 2549970 4629720, "
        "2549970 4629660, 2549940 4629660, 2549939.26359348 4629660, 2549939.26359348 "
        "4629720)), ((2549970 4629780, 2550000 4629780, 2550090 4629780, "
        "2550090 4629790.803233128, 2550092.13765005 4629790.803233128, "
        "2550092.13765005 4629750, 2550090 4629750, 2550090 4629720, "
        "2550092.13765005 4629720, 2550092.13765005 4629660, 2550090 4629660, "
        "2550090 4629637.929176555, 2550060 4629637.929176555, 2550060 4629660, "
        "2550030 4629660, 2550030 4629690, 2550000 4629690, 2550000 4629720, "
        "2549970 4629720, 2549970 4629780)))"
    )

    built_geometry = TileProcessor.polygonize_array(
        pixel_values, pixel_georeferences, datasource
    )

    clipped_built_geometry = TileProcessor.clip_to_tile_extent(built_geometry, tile)

    assert str(clipped_built_geometry) == expected_clipped_geometry