Commit 6aa76f4d authored by Nicolas Garcia Ospina's avatar Nicolas Garcia Ospina
Browse files

Included grid generator scripts

parent de53b79c
Pipeline #25744 passed with stage
in 2 minutes and 29 seconds
......@@ -9,6 +9,8 @@ cache:
- venv/
before_script:
- apt-get update -y
- apt-get install -y python3-pip proj-bin libgeos-dev gdal-bin libgdal-dev
- python3 -V
- pip3 install virtualenv
- virtualenv 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 sys
import geopandas as gpd
import argparse
from babelgrid import Babel
from multiprocessing import Pool
import numpy as np
from shapely.geometry import Polygon
logger = logging.getLogger()
logger.setLevel(logging.INFO)
logger.addHandler(logging.StreamHandler(sys.stdout))
logger.info("Launching quadtreegrid")
parser = argparse.ArgumentParser(
description="""Get all quadkeys from given polygon shapefile"""
)
parser.add_argument("-country", type=str, help="Country name in GADM 3.6")
parser.add_argument(
"-zoom",
nargs="?",
const=1,
type=int,
default=18,
help="Set zoom level of output tiles, default=18",
)
parser.add_argument(
"-max_batches",
nargs="?",
const=1,
type=int,
default=100,
help="""Set the maximum amount of batches for multiprocessing""",
)
parser.add_argument(
"-n_cores",
nargs="?",
const=1,
type=int,
default=1,
help="""Set the maximum amount of jobs to spawn in multiprocessing""",
)
parser.add_argument(
"-output",
type=str,
nargs="?",
const=1,
default="./output_quadkeys.txt",
help="""Set the output filepath and extension as
text or CSV filepath, default=./output_quadkeys.txt""",
)
args = parser.parse_args()
# Get a country polygon by name from the GADM dataset and do small buffer
world = gpd.read_file("./gadm36_0.shp").to_crs("epsg:4326")
country = world[world["NAME_0"] == args.country]
logger.info("Setting a buffer on target geometry")
country.geometry.iloc[:] = country.buffer(0.05).values
# Build the batches from griding the input polygon
sw_lon, sw_lat, ne_lon, ne_lat = country.geometry.values[0].bounds
logger.info("Creating batches")
width = abs(ne_lon - sw_lon) / np.sqrt(args.max_batches)
height = abs(ne_lat - sw_lat) / np.sqrt(args.max_batches)
cols = list(np.arange(int(np.floor(sw_lon)) - 0.5, int(np.ceil(ne_lon)) + 0.5, width))
rows = list(np.arange(int(np.floor(sw_lat)) - 0.5, int(np.ceil(ne_lat)) + 0.5, height))
rows.reverse()
polygons = []
for x in cols:
for y in rows:
polygons.append(
Polygon(
[
(x, y),
(x + width, y),
(x + width, y - height),
(x, y - height),
]
)
)
grid = gpd.GeoDataFrame({"geometry": polygons}, crs="epsg:4326")
gridded_polygon = gpd.overlay(grid, country, how="intersection")
batch_list = gridded_polygon.geometry.tolist()
logger.info("{} to be processed".format(len(batch_list)))
def get_tiles_from_polygon(in_geometry):
"""Wrapper function that extracts quadkeys inside a given polygon.
Args:
in_geometry (shapely.geometry.polygon.Polygon): Polygon to fill with quadkeys.
Returns:
tiles_ids (list): List of quadkeys as strings.
"""
if in_geometry.geom_type == "Polygon":
tiles = Babel("bing").polyfill(in_geometry, resolution=args.zoom)
elif in_geometry.geom_type == "MultiPolygon":
tiles = []
geometries = list(in_geometry)
for geometry in geometries:
tiles.extend(Babel("bing").polyfill(geometry, resolution=args.zoom))
else:
raise IOError("Input is not a geometry.")
tiles_ids = [tile.to_dict()["tile_id"] for tile in tiles]
logger.info("Created list with {} tiles".format(len(tiles)))
return tiles_ids
def main():
logger.info("create pool")
p = Pool(processes=args.n_cores)
logger.info("execute pool")
data = p.map(get_tiles_from_polygon, batch_list)
logger.info("closing pool")
p.close()
p.join()
# Remove possible duplicated quadkeys
flat_list_ids = list(set([item for sublist in data for item in sublist]))
logger.info("Writing {} to {}".format(len(flat_list_ids), args.output))
outfile = open(args.output, "w")
if args.output[-3:] == "csv":
outfile.write("quadkey\n")
for element in flat_list_ids:
outfile.write(element + "\n")
outfile.close()
if __name__ == "__main__":
main()
......@@ -19,18 +19,75 @@
import logging
import sys
import geopandas as gpd
import mercantile
import argparse
# Add a logger printing error, warning, info and debug messages to the screen
logger = logging.getLogger()
logger.setLevel(logging.DEBUG)
logger.addHandler(logging.StreamHandler(sys.stdout))
logger.info("Launching quadtreegrid")
parser = argparse.ArgumentParser(
description="""Get all quadkeys from given bounds
or get the bounds of a vector file"""
)
parser.add_argument(
"-in_bounds",
type=float,
nargs="+",
help="Bounds in order West,South,East,North",
)
parser.add_argument("-in_file", type=str, help="Shapefile with a given geometry")
parser.add_argument(
"-zoom",
nargs="?",
const=1,
type=int,
default=18,
help="Set zoom level of output tiles, default=18",
)
parser.add_argument(
"-output",
type=str,
nargs="?",
const=1,
default="./output_quadkeys.txt",
help="""Set the output filepath and extension as
text or CSV filepath, default=./output_quadkeys.txt""",
)
args = parser.parse_args()
def main():
if args.in_bounds:
sw_lon, sw_lat, ne_lon, ne_lat = args.in_bounds
print(ne_lat)
elif args.in_file:
outline = gpd.read_file(args.in_file).to_crs("epsg:4326")
geometry = outline.envelope.values[0]
sw_lon, sw_lat, ne_lon, ne_lat = geometry.bounds
else:
raise ValueError("No input geometry found")
# Example logging output
logger.info("quadtreegrid has started")
logger.info(
"Zoom level {} quadkeys will be created for the bounding "
"box {} {} {} {}.".format(args.zoom, sw_lon, sw_lat, ne_lon, ne_lat)
)
tiles_generator = mercantile.tiles(sw_lon, sw_lat, ne_lon, ne_lat, args.zoom)
logger.info("Writing to {}".format(args.output))
outfile = open(args.output, "w")
if args.output[-3:] == "csv":
outfile.write("quadkey\n")
for tile in tiles_generator:
outfile.write(mercantile.quadkey(tile) + "\n")
outfile.close()
# Example logging output
logger.info("quadtreegrid is closing")
# Leave the program
sys.exit()
......
......@@ -8,10 +8,17 @@ linters_require = ["black>=20.8b1", "pylint", "flake8"]
setup(
name="quadtreegrid",
version="0.0.1",
description="Quadtree grid creation for OpenBuildingMap",
description="Quadtree grid creation",
keywords="OpenBuildingMap, quadtree, grid",
author="Helmholtz-Zentrum Potsdam Deutsches GeoForschungsZentrum GFZ",
install_requires=["numpy"],
install_requires=[
"pyproj==3.0.0.post1",
"shapely",
"mercantile",
"geopandas",
"numpy",
"babelgrid",
],
extras_require={
"tests": tests_require,
"linters": linters_require,
......
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