Commit a2e55b7b authored by Cecilia Nievas's avatar Cecilia Nievas
Browse files

Rearranged and relocated contents of data_unit_tiles.py

parent 8b42395e
Pipeline #36320 passed with stage
in 1 minute and 43 seconds
......@@ -18,13 +18,11 @@
import logging
import numpy
import pandas
import geopandas
import pyproj
import mercantile
import shapely
from copy import deepcopy
from gdeimporter.tools.database import Database
from gdeimporter.tools.spatial import SpatialTools
logger = logging.getLogger()
......@@ -66,7 +64,7 @@ class DataUnitTilesHelper:
"""
# Retrieve all quadtiles associated with the bounding box of 'in_geometry'
quadtiles_bbox = DataUnitTilesHelper.get_quadtiles_in_bbox_of_polygon(in_geometry)
quadtiles_bbox = SpatialTools.get_quadtiles_in_bbox_of_polygon(in_geometry)
# Create a GeoPandas GeoDataFrame with the original geometry
in_geometry_df = geopandas.GeoDataFrame({"geometry": [in_geometry]})
......@@ -187,7 +185,7 @@ class DataUnitTilesHelper:
all_fractions = numpy.zeros([data_unit_tiles.shape[0]])
# Tiles that are fully contained in the original geometry 'data_unit_geometry'
all_areas[data_unit_tiles.contained] = DataUnitTilesHelper.get_area_on_sphere(
all_areas[data_unit_tiles.contained] = SpatialTools.get_area_on_sphere(
data_unit_tiles[data_unit_tiles.contained]["lon_w"].values,
data_unit_tiles[data_unit_tiles.contained]["lat_s"].values,
data_unit_tiles[data_unit_tiles.contained]["lon_e"].values,
......@@ -240,146 +238,6 @@ class DataUnitTilesHelper:
return data_unit_tiles
@staticmethod
def get_quadtiles_in_bbox_of_polygon(in_geometry, zoomlevel=18):
"""This function retrieves the quadtiles of zoom level 'zoomlevel' that fully contain
the bounding box(es) of 'in_geometry'. If 'in_geometry' is a MultiPolygon, the bounding
boxes are defined for each Polygon in the MultiPolygon, not the MultiPolygon as a whole.
The quadtiles are defined in EPSG:3857 projection, but their geometry is given in the
WGS84 coordinate system (ESPG:4326).
Args:
in_geometry (Shapely Polygon or MultiPolygon):
Geometry for which the associated quadtiles will be retrieved, defined in
ESPG:4326.
zoomlevel (int):
Zoom level of the quadtiles to be retrieved.
Returns:
grid (GeoPandas GeoDataFrame):
GeoDataFrame containing the definition of the quadtiles in terms of:
quadkey (str):
Quadkey of the level 18 quadtile.
lon_w, lat_s, lon_e, lat_n (float)
Coordinates of the corners of the quadtile.
geometry (Shapely Polygon):
Geometry of the quadtile.
"""
if in_geometry.geom_type == "Polygon":
geometries = [in_geometry]
elif in_geometry.geom_type == "MultiPolygon":
geometries = list(in_geometry)
else:
raise OSError("get_quadtiles_in_bbox_of_polygon: argument is not a valid geometry")
tiles_quadkeys = []
tiles_lons_w = []
tiles_lats_s = []
tiles_lons_e = []
tiles_lats_n = []
for geometry in geometries:
if DataUnitTilesHelper.crosses_antimeridian(geometry):
raise OSError(
"get_quadtiles_in_bbox_of_polygon: geometry crosses the antimeridian"
)
lon_w, lat_s, lon_e, lat_n = geometry.bounds
if lat_s < -85.0511 or lat_n > 85.0511:
raise OSError(
"get_quadtiles_in_bbox_of_polygon: "
"quadtiles cannot be defined for latitudes outside [-85.0511, +85.0511]"
)
tiles = list(mercantile.tiles(lon_w, lat_s, lon_e, lat_n, zoomlevel))
quadkeys = list([mercantile.quadkey(tile) for tile in tiles])
lons_w = list([mercantile.bounds(tile).west for tile in tiles])
lats_s = list([mercantile.bounds(tile).south for tile in tiles])
lons_e = list([mercantile.bounds(tile).east for tile in tiles])
lats_n = list([mercantile.bounds(tile).north for tile in tiles])
tiles_quadkeys.extend(quadkeys)
tiles_lons_w.extend(lons_w)
tiles_lats_s.extend(lats_s)
tiles_lons_e.extend(lons_e)
tiles_lats_n.extend(lats_n)
# Create a Pandas DataFrame
aux_grid = pandas.DataFrame(
{
"quadkey": tiles_quadkeys,
"lon_w": tiles_lons_w,
"lat_s": tiles_lats_s,
"lon_e": tiles_lons_e,
"lat_n": tiles_lats_n,
}
)
aux_grid = aux_grid.drop_duplicates(ignore_index=True) # Remove potential duplicates
grid = DataUnitTilesHelper.create_geometry_of_quadtiles(aux_grid)
grid = grid.set_index("quadkey", drop=False)
return grid
@staticmethod
def create_geometry_of_quadtiles(in_grid):
"""This function generates the geometry of the tiles defined by their corners in
'in_grid' and outputs the latter (in ESPG:4326) together with all the contents of
'in_grid' as a GeoPandas GeoDataFrame.
Args:
in_grid (Pandas DataFrame):
DataFrame with at least the following columns:
lon_w, lat_s, lon_e, lat_n (float):
Coordinates of the corners of the quadtile (in degrees).
Returns:
out_grid (GeoPandas GeoDataFrame):
GeoDataFrame with the same contents as in_grid plus the geometry of the
quadtiles in ESPG:4326.
"""
geoms = []
for i in range(in_grid.shape[0]):
lon_w = in_grid["lon_w"].values[i]
lat_s = in_grid["lat_s"].values[i]
lon_e = in_grid["lon_e"].values[i]
lat_n = in_grid["lat_n"].values[i]
geoms.append(
shapely.geometry.Polygon(
[(lon_w, lat_s), (lon_e, lat_s), (lon_e, lat_n), (lon_w, lat_n)]
)
)
geometries = geopandas.GeoSeries(geoms)
out_grid = geopandas.GeoDataFrame(in_grid, geometry=geometries)
out_grid.crs = pyproj.CRS("epsg:4326")
return out_grid
@staticmethod
def crosses_antimeridian(in_geometry):
"""This function checks whether 'in_geometry' crosses the antimeridian (i.e. longitude
equal to -180 or +180.
Args:
in_geometry (Shapely Polygon or MultiPolygon):
Geometry for which the function will assess if the antimeridian is crossed,
defined in ESPG:4326.
Returns (bool):
True if the antimeridian is crossed, False if the antimeridian is not crossed.
"""
longitude_width = abs(in_geometry.bounds[2] - in_geometry.bounds[0])
if longitude_width < 180.0:
return False
else:
return True
@staticmethod
def determine_if_areas_are_the_same(data_unit_tiles, tiles):
"""This function compares the areas of the geometries of 'data_unit_tiles' and 'tiles'
......@@ -418,51 +276,12 @@ class DataUnitTilesHelper:
return are_the_same
@staticmethod
def get_area_on_sphere(lon_w, lat_s, lon_e, lat_n, earth_sphere_radius=6371000.0):
"""This function returns the area (m2) of a rectangle on Earth calculated assuming the
Earth is a sphere with radius 'earth_sphere_radius'. The rectangle is defined by its
four bounding coordinates 'lon_w', 'lat_s', 'lon_e' and 'lat_n', all given in degrees.
When used to calculate the area of zoom level 18 quadtiles, this approximate function
yields area values that differ in less than 1% with respect to those calculated using
the Albers equal area projection.
References:
Kelly, K., Šavrič, B. (2021). Area and volume computation of longitude–latitude
grids and three-dimensional meshes. Transactions in GIS, 25(1):6–24.
Args:
lon_w (float or array of floats):
West bound of the tile(s), in degrees.
lat_s (float or array of floats):
South bound of the tile(s), in degrees.
lon_e (float or array of floats):
East bound of the tile(s), in degrees.
lat_n (float or array of floats):
North bound of the tile(s), in degrees.
earth_sphere_radius (float):
Radius of the Earth assuming a sphere, in meters. Default: 6,371,000 m.
Returns:
area (float or array of floats):
Area of the rectangles defined by the input coordinates, in square metres.
"""
area = (
(numpy.pi / 180.0)
* numpy.power(earth_sphere_radius, 2)
* numpy.abs(numpy.sin(numpy.deg2rad(lat_n)) - numpy.sin(numpy.deg2rad(lat_s)))
* numpy.abs(lon_e - lon_w)
)
return area
@staticmethod
def get_areas_and_ratios_arbitrary_shapes(geometries_1, geometries_2):
"""This function uses get_area_albers_equal_area() to calculate the area of each polygon
of 'geometries_1' and 'geometries_2' and calculates the ratio of the former with respect
to the latter. Both 'geometries_1' and 'geometries_2' need to have the same length.
"""This function uses SpatialTools.get_area_albers_equal_area() to calculate the area of
each polygon of 'geometries_1' and 'geometries_2' and calculates the ratio of the former
with respect to the latter. Both 'geometries_1' and 'geometries_2' need to have the same
length.
Args:
geometries_1, geometries_2 (arrays of Shapely Polygons):
......@@ -480,42 +299,13 @@ class DataUnitTilesHelper:
areas_2 = numpy.zeros([len(geometries_1)])
for i in range(len(geometries_1)):
areas_1[i] = DataUnitTilesHelper.get_area_albers_equal_area(geometries_1[i])
areas_2[i] = DataUnitTilesHelper.get_area_albers_equal_area(geometries_2[i])
areas_1[i] = SpatialTools.get_area_albers_equal_area(geometries_1[i])
areas_2[i] = SpatialTools.get_area_albers_equal_area(geometries_2[i])
ratios = areas_1 / areas_2
return areas_1, areas_2, ratios
@staticmethod
def get_area_albers_equal_area(input_polygon):
"""This function returns the area (m2) of a polygon on Earth calculated using the Albers
equal area projection. The polygon's geometry is assumed to be defined in ESPG:4326.
Args:
in_geometry (Shapely polygon):
Polygon defined in ESPG:4326.
Returns:
projected_area: Area of the polygon in m2 (in Albers equal area projection).
"""
# Get the bounding box of the polygon
lon_w = input_polygon.bounds[0]
lat_s = input_polygon.bounds[1]
lon_e = input_polygon.bounds[2]
lat_n = input_polygon.bounds[3]
# Use the resulting coordinates to carry out the transformation
project_albers = pyproj.Proj(
"+proj=aea +lat_1={} +lat_2={} +lat_0={} +lon_0={}".format(
lat_s, lat_n, (lat_s + lat_n) / 2.0, (lon_w + lon_e) / 2.0
)
)
geometry = shapely.ops.transform(project_albers, input_polygon)
return geometry.area
@staticmethod
def retrieve_built_up_area(quadkeys, db_built_up_config, db_table):
"""This function retrieves the built-up area associated with the input quadkeys from
......
......@@ -20,7 +20,7 @@ import logging
import iso3166
from multiprocessing import Pool
from functools import partial
from gdeimporter.tools.data_unit_tiles import DataUnitTilesHelper
from gdeimporter.data_unit_tiles import DataUnitTilesHelper
from gdeimporter.tools.database import Database
......
#!/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 numpy
import pandas
import geopandas
import pyproj
import mercantile
import shapely
logger = logging.getLogger()
class SpatialTools:
"""This class contains methods that carry out spatial operations."""
@staticmethod
def get_quadtiles_in_bbox_of_polygon(in_geometry, zoomlevel=18):
"""This function retrieves the quadtiles of zoom level 'zoomlevel' that fully contain
the bounding box(es) of 'in_geometry'. If 'in_geometry' is a MultiPolygon, the bounding
boxes are defined for each Polygon in the MultiPolygon, not the MultiPolygon as a whole.
The quadtiles are defined in EPSG:3857 projection, but their geometry is given in the
WGS84 coordinate system (ESPG:4326).
Args:
in_geometry (Shapely Polygon or MultiPolygon):
Geometry for which the associated quadtiles will be retrieved, defined in
ESPG:4326.
zoomlevel (int):
Zoom level of the quadtiles to be retrieved.
Returns:
grid (GeoPandas GeoDataFrame):
GeoDataFrame containing the definition of the quadtiles in terms of:
quadkey (str):
Quadkey of the level 18 quadtile.
lon_w, lat_s, lon_e, lat_n (float)
Coordinates of the corners of the quadtile.
geometry (Shapely Polygon):
Geometry of the quadtile.
"""
if in_geometry.geom_type == "Polygon":
geometries = [in_geometry]
elif in_geometry.geom_type == "MultiPolygon":
geometries = list(in_geometry)
else:
raise OSError("get_quadtiles_in_bbox_of_polygon: argument is not a valid geometry")
tiles_quadkeys = []
tiles_lons_w = []
tiles_lats_s = []
tiles_lons_e = []
tiles_lats_n = []
for geometry in geometries:
if SpatialTools.crosses_antimeridian(geometry):
raise OSError(
"get_quadtiles_in_bbox_of_polygon: geometry crosses the antimeridian"
)
lon_w, lat_s, lon_e, lat_n = geometry.bounds
if lat_s < -85.0511 or lat_n > 85.0511:
raise OSError(
"get_quadtiles_in_bbox_of_polygon: "
"quadtiles cannot be defined for latitudes outside [-85.0511, +85.0511]"
)
tiles = list(mercantile.tiles(lon_w, lat_s, lon_e, lat_n, zoomlevel))
quadkeys = list([mercantile.quadkey(tile) for tile in tiles])
lons_w = list([mercantile.bounds(tile).west for tile in tiles])
lats_s = list([mercantile.bounds(tile).south for tile in tiles])
lons_e = list([mercantile.bounds(tile).east for tile in tiles])
lats_n = list([mercantile.bounds(tile).north for tile in tiles])
tiles_quadkeys.extend(quadkeys)
tiles_lons_w.extend(lons_w)
tiles_lats_s.extend(lats_s)
tiles_lons_e.extend(lons_e)
tiles_lats_n.extend(lats_n)
# Create a Pandas DataFrame
aux_grid = pandas.DataFrame(
{
"quadkey": tiles_quadkeys,
"lon_w": tiles_lons_w,
"lat_s": tiles_lats_s,
"lon_e": tiles_lons_e,
"lat_n": tiles_lats_n,
}
)
aux_grid = aux_grid.drop_duplicates(ignore_index=True) # Remove potential duplicates
grid = SpatialTools.create_geometry_of_quadtiles(aux_grid)
grid = grid.set_index("quadkey", drop=False)
return grid
@staticmethod
def create_geometry_of_quadtiles(in_grid):
"""This function generates the geometry of the tiles defined by their corners in
'in_grid' and outputs the latter (in ESPG:4326) together with all the contents of
'in_grid' as a GeoPandas GeoDataFrame.
Args:
in_grid (Pandas DataFrame):
DataFrame with at least the following columns:
lon_w, lat_s, lon_e, lat_n (float):
Coordinates of the corners of the quadtile (in degrees).
Returns:
out_grid (GeoPandas GeoDataFrame):
GeoDataFrame with the same contents as in_grid plus the geometry of the
quadtiles in ESPG:4326.
"""
geoms = []
for i in range(in_grid.shape[0]):
lon_w = in_grid["lon_w"].values[i]
lat_s = in_grid["lat_s"].values[i]
lon_e = in_grid["lon_e"].values[i]
lat_n = in_grid["lat_n"].values[i]
geoms.append(
shapely.geometry.Polygon(
[(lon_w, lat_s), (lon_e, lat_s), (lon_e, lat_n), (lon_w, lat_n)]
)
)
geometries = geopandas.GeoSeries(geoms)
out_grid = geopandas.GeoDataFrame(in_grid, geometry=geometries)
out_grid.crs = pyproj.CRS("epsg:4326")
return out_grid
@staticmethod
def crosses_antimeridian(in_geometry):
"""This function checks whether 'in_geometry' crosses the antimeridian (i.e. longitude
equal to -180 or +180.
Args:
in_geometry (Shapely Polygon or MultiPolygon):
Geometry for which the function will assess if the antimeridian is crossed,
defined in ESPG:4326.
Returns (bool):
True if the antimeridian is crossed, False if the antimeridian is not crossed.
"""
longitude_width = abs(in_geometry.bounds[2] - in_geometry.bounds[0])
if longitude_width < 180.0:
return False
else:
return True
@staticmethod
def get_area_on_sphere(lon_w, lat_s, lon_e, lat_n, earth_sphere_radius=6371000.0):
"""This function returns the area (m2) of a rectangle on Earth calculated assuming the
Earth is a sphere with radius 'earth_sphere_radius'. The rectangle is defined by its
four bounding coordinates 'lon_w', 'lat_s', 'lon_e' and 'lat_n', all given in degrees.
When used to calculate the area of zoom level 18 quadtiles, this approximate function
yields area values that differ in less than 1% with respect to those calculated using
the Albers equal area projection.
References:
Kelly, K., Šavrič, B. (2021). Area and volume computation of longitude–latitude
grids and three-dimensional meshes. Transactions in GIS, 25(1):6–24.
Args:
lon_w (float or array of floats):
West bound of the tile(s), in degrees.
lat_s (float or array of floats):
South bound of the tile(s), in degrees.
lon_e (float or array of floats):
East bound of the tile(s), in degrees.
lat_n (float or array of floats):
North bound of the tile(s), in degrees.
earth_sphere_radius (float):
Radius of the Earth assuming a sphere, in meters. Default: 6,371,000 m.
Returns:
area (float or array of floats):
Area of the rectangles defined by the input coordinates, in square metres.
"""
area = (
(numpy.pi / 180.0)
* numpy.power(earth_sphere_radius, 2)
* numpy.abs(numpy.sin(numpy.deg2rad(lat_n)) - numpy.sin(numpy.deg2rad(lat_s)))
* numpy.abs(lon_e - lon_w)
)
return area
@staticmethod
def get_area_albers_equal_area(input_polygon):
"""This function returns the area (m2) of a polygon on Earth calculated using the Albers
equal area projection. The polygon's geometry is assumed to be defined in ESPG:4326.
Args:
in_geometry (Shapely polygon):
Polygon defined in ESPG:4326.
Returns:
projected_area: Area of the polygon in m2 (in Albers equal area projection).
"""
# Get the bounding box of the polygon
lon_w = input_polygon.bounds[0]
lat_s = input_polygon.bounds[1]
lon_e = input_polygon.bounds[2]
lat_n = input_polygon.bounds[3]
# Use the resulting coordinates to carry out the transformation
project_albers = pyproj.Proj(
"+proj=aea +lat_1={} +lat_2={} +lat_0={} +lon_0={}".format(
lat_s, lat_n, (lat_s + lat_n) / 2.0, (lon_w + lon_e) / 2.0
)
)
geometry = shapely.ops.transform(project_albers, input_polygon)
return geometry.area
......@@ -22,10 +22,9 @@ import numpy
import geopandas
import pandas
import shapely
import pytest
from copy import deepcopy
from pyproj import CRS
from gdeimporter.tools.data_unit_tiles import DataUnitTilesHelper
from gdeimporter.data_unit_tiles import DataUnitTilesHelper
from gdeimporter.configuration import Configuration
from gdeimporter.tools.database import Database
from gdeimporter.dataunit import DataUnit
......@@ -33,179 +32,6 @@ from gdeimporter.dataunit import DataUnit
logger = logging.getLogger()
def test_get_quadtiles_in_bbox_of_polygon():
# Read shapefile with input polygons to test
input_polygons = geopandas.read_file(
os.path.join(
os.path.dirname(__file__),
"data",
"polygons_for_testing_data_unit_tiles.shp",
)
)
# Define expected output:
expected_quadtiles = {}
# Test case 1: Normal polygon in Luxembourg
expected_quadtiles[1] = [
"120203202122231123",
"120203202122231301",
"120203202122231303",
"120203202122231132",
"120203202122231310",
"120203202122231312",
"120203202122231133",
"120203202122231311",
"120203202122231313",
"120203202122320022",
"120203202122320200",
"120203202122320202",
]
# Test case 2: polygon smaller than a zoom level 18 tile in Luxembourg
expected_quadtiles[2] = ["120203202122320210"]
# Test case 3: polygon with west bound = -180.0, northern hemisphere
expected_quadtiles[3] = [
"020202202022220200",
"020202202022220202",
"020202202022220220",
"020202202022220201",
"020202202022220203",
"020202202022220221",
]
# Test case 4: polygon with west bound = -180.0, southern hemisphere
expected_quadtiles[4] = [
"202020020200002002",
"202020020200002020",
"202020020200002022",
"202020020200002003",
"202020020200002021",
"202020020200002023",
]
# Test case 5: polygon with east bound = +180.0, southern hemisphere
expected_quadtiles[5] = [
"313131131311113112",
"313131131311113130",
"313131131311113132",
"313131131311113113",
"313131131311113131",
"313131131311113133",
]
# Test case 6: polygon with east bound = +180.0, northern hemisphere
expected_quadtiles[6] = [
"131313313133331310",
"131313313133331312",
"131313313133331330",
"131313313133331311",
"131313313133331313",
"131313313133331331",
]
# Test case 7: polygon crossing Greenwich, northern hemisphere
expected_quadtiles[7] = [
"031313131131313313",