Commit 8404c6ab authored by Cecilia Nievas's avatar Cecilia Nievas
Browse files

Re-wrote GDE_TOOLS_world_grid.py to work with Quadtiles level 18

parent 5999f08a
"""
Copyright (C) 2020
Helmholtz-Zentrum Potsdam Deutsches GeoForschungsZentrum GFZ
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
......@@ -23,17 +23,11 @@ Section 2.6: Seismic Hazard and Risk Dynamics
GDE_TOOLS_world_grid
====================
THIS CODE ALLOWS TO GENERATE THE WORLD GRID THAT WE
USE FOR OUR DYNAMIC EXPOSURE MODEL
HARD-CODED ASSUMPTIONS ARE:
- The grid spacing is 10 arc-seconds.
- The grid runs from -180.0 through +180.0 in longitude
- The grid runs from -90.0 through +90.0 in latitude
- The top-left-most cell (NW) is cell number 0
- The cell id increases from this first cell to the east, by "row"
- At the end of each row, the cell id "jumps" to the first (westmost)
cell of the next row.
This code allows to generate the world grid that we use for our dynamic exposure model.
It is a zoom level 18 Quadkey-based grid.
We use mercantile (https://pypi.org/project/mercantile/1.2.1/) to create it.
The zoom level (18) is a default in the functions in this file.
"""
import numpy as np
......@@ -42,132 +36,151 @@ import pandas as pd
import shapely
import pyproj
from functools import partial
import mercantile
def adjust_lon(coord_lon):
"""This function adjusts the longitude so that it falls in the range [-180.0, +180.0].
Args:
coord_lon (float): Longitude.
Returns:
coord_lon (float): Longitude in the range: [-180.0,+180.0].
"""
This function adjusts the longitude so that get_cell_id gives the correct cell in the extremes of the grid.
The adjustment is based on the assumption of a 10 arcsec grid.
"""
lon_lowest= -180.0+(5.0/60.0)/60.0
lon_highest= 180.0-(5.0/60.0)/60.0
if coord_lon<lon_lowest:
return lon_lowest
elif coord_lon>lon_highest:
return lon_highest
if coord_lon<0.0:
return max(coord_lon, -180.0)
else:
return coord_lon
return min(coord_lon, 180.0)
def adjust_lat(coord_lat):
""" This function adjusts the latitude so that it falls in the range [-85.0511,+85.0511].
Args:
coord_lat (float): Latitude.
Returns:
coord_lat (float): Latitude in the range: [-85.0511,+85.0511].
"""
This function adjusts the latitude so that get_cell_id gives the correct cell in the extremes of the grid.
The adjustment is based on the assumption of a 10 arcsec grid.
"""
lat_lowest= -90.0+(5.0/60.0)/60.0
lat_highest= 90.0-(5.0/60.0)/60.0
if coord_lat<lat_lowest:
return lat_lowest
elif coord_lat>lat_highest:
return lat_highest
if coord_lat<0.0:
return max(coord_lat, -85.0511)
else:
return coord_lat
return min(coord_lat, 85.0511)
def get_cell_id(coord_lon, coord_lat):
"""
coord_lon= [-180.0,+180.0]
coord_lat= [-90.0,+90.0]
Given any pair of (lon,lat) coordinates, this function returns the ID of the cell
where this point belongs.
This function needs the functions adjust_lon and adjust_lat to make sure that
it gives the correct cell in the extremes of the grid.
def get_cell_id(coord_lon, coord_lat, zoomlevel=18):
"""Given any pair of (lon,lat) coordinates, this function returns the ID of the Quadtile
where this point belongs, at the zoom level indicated by zoomlevel.
If coord_lon and/or coord_lat are outside the range for which Quadtiles are defined (see
below), the mercantile library "brings" them to the closes extreme of the range (e.g.
lon=180.1 is assumed to be lon=180.0.
Args:
coord_lon (float): Longitude in the range: [-180.0,+180.0].
coord_lat (float): Latitude in the range: [-85.0511,+85.0511].
zoomlevel (int): Zoom level of the tile (>=1).
Returns:
cell_id (str): ID of the Quadtile of zoon level zoomlevel where the point defined by
(coord_lon, coord_lat) lies.
"""
lon_edges= np.linspace(-180.0,180.0,num=129601)
lat_edges_inv= np.linspace(-90.0,90.0,num=64801)
cell_x= np.searchsorted(lon_edges,adjust_lon(coord_lon)) - 1 # position of the cell within the longitude vector
cell_y= np.searchsorted(lat_edges_inv,adjust_lat(coord_lat)) - 1 # position of the cell within the latitude vector
cell_id= cell_x + 8397950400 - 129600 * cell_y
if zoomlevel < 1:
print("ERROR in get_cells_in_bbox: zoomlevel needs to be equal to or larger than 1")
return []
cell_id = mercantile.quadkey(mercantile.tile(coord_lon,
coord_lat,
zoomlevel))
return cell_id
def get_corner_cells_of_bbox(sw_lon, sw_lat, ne_lon, ne_lat):
"""
This function returns the IDs of the four cells that make up the
corner of the bounding box defined by the input coordinates.
sw_lon= longitude of the South-West corner
sw_lat= latitude of the South-West corner
ne_lon= longitude of the North-East corner
ne_lat= latitude of the North-East corner
The output order is: South-West corner, South-East corner, North-East corner, North-West corner
(i.e. counter-clockwise, starting from South-west corner).
"""
cell_sw_id= get_cell_id(sw_lon, sw_lat)
cell_ne_id= get_cell_id(ne_lon, ne_lat)
cell_nw_id= get_cell_id(sw_lon, ne_lat)
cell_se_id= get_cell_id(ne_lon, sw_lat)
return cell_sw_id, cell_se_id, cell_ne_id, cell_nw_id
def get_cells_in_bbox(sw_lon, sw_lat, ne_lon, ne_lat):
"""
This function returns the IDs of all the cells contained within the bounding box
defined by the input coordinates.
sw_lon= longitude of the South-West corner
sw_lat= latitude of the South-West corner
ne_lon= longitude of the North-East corner
ne_lat= latitude of the North-East corner
The output list starts with the North-West corner and moves first to the East, to then
jump to the row to the South, and then goes again from West to East.
def get_cells_in_bbox(sw_lon, sw_lat, ne_lon, ne_lat, zoomlevel=18):
"""This function returns the IDs of all the cells contained within the bounding box defined
by the input coordinates.
Args:
sw_lon (float): Longitude of the south-west corner.
sw_lat (float): Latitude of the south-west corner.
ne_lon (float): Longitude of the north-east corner.
ne_lat (float): Latitude of the north-east corner.
zoomlevel (int): Zoom level of the tile (>=1).
Returns:
cells_ids (list of strings)
"""
cell_sw_id, cell_se_id, cell_ne_id, cell_nw_id= get_corner_cells_of_bbox(sw_lon, sw_lat, ne_lon, ne_lat)
ew_width= cell_ne_id - (cell_nw_id - 1) # number of cells in the width of the bbox
ns_height= int((cell_sw_id - cell_nw_id) / 129600 + 1) # number of cells in the height of the bbox
cells_ids= []
count= 0
for i in range(0,ns_height):
for j in range(0,ew_width):
if ew_width*ns_height>50000:
if count%50000==0:
print(' Working on cell '+str(count+1)+' of '+str(ew_width*ns_height))
count= count + 1
cells_ids.append(cell_nw_id+i*129600+j)
if sw_lon > ne_lon:
print("ERROR in get_cells_in_bbox: sw_lon should be smaller than ne_lon")
return []
if sw_lat > ne_lat:
print("ERROR in get_cells_in_bbox: sw_lat should be smaller than ne_lat")
return []
if zoomlevel < 1:
print("ERROR in get_cells_in_bbox: zoomlevel needs to be equal to or larger than 1")
return []
all_tiles = mercantile.tiles(adjust_lon(sw_lon),
adjust_lat(sw_lat),
adjust_lon(ne_lon),
adjust_lat(ne_lat),
zoomlevel)
cells_ids = []
for tile in all_tiles:
cells_ids.append(mercantile.quadkey(tile))
return cells_ids
def get_coordinates_from_cell_id(cell_id):
"""Given a cell ID, it returns the coordinates that define the four corners of the tile.
Args:
cell_id (str): Quadkey of the tile. The number of digits in the string implicitly
indicate the associated zoom level.
Returns:
lon_w (float): West longitude of the tile.
lat_s (float): South latitude of the tile.
lon_e (float): East longitude of the tile.
lat_n (float): North latitude of the tile.
"""
Given a cell ID, it returns the coordinates that define the four corners of the cell.
"""
if cell_id<8398080000 and cell_id>-1:
lon_edges= np.linspace(-180.0,180.0,num=129601)
lat_edges_inv= np.linspace(-90.0,90.0,num=64801)
cell_y= -(cell_id // 129600) +64799
cell_x= cell_id - 8397950400 + 129600 * cell_y
lon_w= lon_edges[cell_x]
lon_e= lon_edges[cell_x+1]
lat_s= lat_edges_inv[cell_y]
lat_n= lat_edges_inv[cell_y+1]
else:
print('ERROR!! CELL ID '+str(cell_id)+' DOES NOT EXIST!!!')
lon_w= -999.9
lon_e= -999.9
lat_s= -999.9
lat_n= -999.9
tile_bounds = mercantile.bounds(mercantile.quadkey_to_tile(cell_id))
lon_w = tile_bounds.west
lon_e = tile_bounds.east
lat_s = tile_bounds.south
lat_n = tile_bounds.north
return lon_w, lat_s, lon_e, lat_n
def get_cells_geometry_as_geopandas(cell_id_list, filename=None):
"""This function returns a GeoPandas DataFrame of all tiles in cell_id_list.
If a filename is given, it will export it as a GeoJSON.
Args:
cell_id_list (list of strings): List of tile IDs.
Returns:
out_gdf (GeoDataFrame): GeoDataFrame with the tile IDs and another with their
geometries.
"""
This function returns a GeoPandas DataFrame of all cells in cell_id_list.
If a filename is given, it will export it as a ShapeFile.
cell_id_list= list of cell IDs
"""
geom_list= []
geom_list = []
for cell_id in cell_id_list:
lon_w, lat_s, lon_e, lat_n= get_coordinates_from_cell_id(cell_id)
geom_list.append(shapely.geometry.Polygon([(lon_w,lat_s), (lon_e,lat_s), (lon_e,lat_n), (lon_w,lat_n)]))
geom_list.append(shapely.geometry.Polygon([(lon_w, lat_s),
(lon_e, lat_s),
(lon_e, lat_n),
(lon_w, lat_n)]))
geometry = gpd.GeoSeries(geom_list)
data= {'cell_id': cell_id_list,
'geometry': geometry}
......@@ -176,135 +189,208 @@ def get_cells_geometry_as_geopandas(cell_id_list, filename=None):
if filename:
print(' Saving to GeoJSON...')
out_gdf.to_file(filename,driver='GeoJSON')
return out_gdf
def get_geometry_WKT_of_cell(cell_id):
"""This function returns the Well-Known Text of the geometry of the tile with cell_id.
"""
lon_w, lat_s, lon_e, lat_n= get_coordinates_from_cell_id(cell_id)
out_wkt= (shapely.geometry.Polygon([(lon_w,lat_s), (lon_e,lat_s), (lon_e,lat_n), (lon_w,lat_n)])).wkt
out_wkt= (shapely.geometry.Polygon([(lon_w, lat_s),
(lon_e, lat_s),
(lon_e, lat_n),
(lon_w, lat_n)])).wkt
return out_wkt
def get_cells_geometry_as_pandas_with_WKT(cell_id_list):
""" This function returns a Pandas DataFrame of all tiles in cell_id_list.
The geometry is returned as Well-Known Text.
Args:
cell_id_list (list of strings): List of tile IDs.
Returns:
out_df (DataFrame): DataFrame with a column for the tile IDs and another for their
geometries in Well-Known Text format.
"""
This function returns a Pandas DataFrame of all cells in cell_id_list.
If a filename is given, it will export it as a ShapeFile.
The geometry is given as Well Known Text.
cell_id_list= list of cell IDs
"""
geom_list= []
geom_list = []
for cell_id in cell_id_list:
lon_w, lat_s, lon_e, lat_n= get_coordinates_from_cell_id(cell_id)
geom_list.append((shapely.geometry.Polygon([(lon_w,lat_s), (lon_e,lat_s), (lon_e,lat_n), (lon_w,lat_n)])).wkt)
data= {'cell_id': cell_id_list,
geom_list.append((shapely.geometry.Polygon([(lon_w, lat_s),
(lon_e, lat_s),
(lon_e, lat_n),
(lon_w, lat_n)])).wkt)
data = {'cell_id': cell_id_list,
'geometry': geom_list}
out_df= pd.DataFrame(data, columns = ['cell_id','geometry'])
out_df = pd.DataFrame(data, columns = ['cell_id','geometry'])
return out_df
def get_cells_in_bbox_as_geopandas(sw_lon, sw_lat, ne_lon, ne_lat, zoomlevel=18, filename=None):
"""This function returns a GeoPandas DataFrame of all tiles contained within the bounding
box defined by the input coordinates. If a filename is given, it will export it as a
ShapeFile.
Args:
sw_lon (float): Longitude of the south-west corner.
sw_lat (float): Latitude of the south-west corner.
ne_lon (float): Longitude of the north-east corner.
ne_lat (float): Latitude of the north-east corner.
zoomlevel (int): Zoom level of the tile (>=1).
filename (str): Filename for exporting the GeoDataFrame as a Shapefile.
def get_cells_in_bbox_as_geopandas(sw_lon, sw_lat, ne_lon, ne_lat, filename=None):
"""
This function returns a GeoPandas DataFrame of all cells contained within the bounding box
defined by the input coordinates. If a filename is given, it will export it as a ShapeFile.
sw_lon= longitude of the South-West corner
sw_lat= latitude of the South-West corner
ne_lon= longitude of the North-East corner
ne_lat= latitude of the North-East corner
The output starts with the North-West corner and moves first to the East, to then
jump to the row to the South, and then goes again from West to East.
Returns:
out_gdf (GeoDataFrame): GeoDataFrame with the tile IDs and another with their
geometries.
"""
print(' Gathering cell IDs...')
all_cells_ids= get_cells_in_bbox(sw_lon, sw_lat, ne_lon, ne_lat)
geom_list= []
ids_list= []
print(' Gathering geometry of cells...')
for i, id_i in enumerate(all_cells_ids):
if len(all_cells_ids)>50000:
if i%50000==0:
print(' Working on cell '+str(i+1)+' of '+str(len(all_cells_ids)))
lon_w, lat_s, lon_e, lat_n= get_coordinates_from_cell_id(id_i)
geom_list.append(shapely.geometry.Polygon([(lon_w,lat_s), (lon_e,lat_s), (lon_e,lat_n), (lon_w,lat_n)]))
ids_list.append(str(id_i))
if sw_lon > ne_lon:
print("ERROR in get_cells_in_bbox: sw_lon should be smaller than ne_lon")
return []
if sw_lat > ne_lat:
print("ERROR in get_cells_in_bbox: sw_lat should be smaller than ne_lat")
return []
if zoomlevel < 1:
print("ERROR in get_cells_in_bbox: zoomlevel needs to be equal to or larger than 1")
return []
all_tiles = mercantile.tiles(adjust_lon(sw_lon),
adjust_lat(sw_lat),
adjust_lon(ne_lon),
adjust_lat(ne_lat),
zoomlevel)
cells_ids = []
geom_list = []
print(' Gathering IDs and geometries of tiles...')
for tile in all_tiles:
cells_ids.append(mercantile.quadkey(tile))
tile_bounds = mercantile.bounds(tile)
geom_list.append(shapely.geometry.Polygon([(tile_bounds.west, tile_bounds.south),
(tile_bounds.east, tile_bounds.south),
(tile_bounds.east, tile_bounds.north ),
(tile_bounds.west, tile_bounds.north )]))
print(' Creating GeoPandas DataFrame...')
geometry = gpd.GeoSeries(geom_list)
data= {'cell_id': ids_list,
data= {'cell_id': cells_ids,
'geometry': geometry}
out_gdf= gpd.GeoDataFrame(data, columns = ['cell_id','geometry'])
out_gdf.crs = {'init' :'epsg:4326'}
if filename:
print(' Saving to ShapeFile...')
out_gdf.to_file(filename,driver='ESRI Shapefile')
return out_gdf
def get_cells_in_bbox_as_pandas_with_WKT(sw_lon, sw_lat, ne_lon, ne_lat):
"""
This function returns a Pandas DataFrame of all cells contained within the bounding box
defined by the input coordinates. The geometry is given as Well Known Text.
sw_lon= longitude of the South-West corner
sw_lat= latitude of the South-West corner
ne_lon= longitude of the North-East corner
ne_lat= latitude of the North-East corner
def get_cells_in_bbox_as_pandas_with_WKT(sw_lon, sw_lat, ne_lon, ne_lat, zoomlevel=18):
"""This function returns a Pandas DataFrame of all tiles contained within the bounding
box defined by the input coordinates. The geometry is given as Well Known Text.
Args:
sw_lon (float): Longitude of the south-west corner.
sw_lat (float): Latitude of the south-west corner.
ne_lon (float): Longitude of the north-east corner.
ne_lat (float): Latitude of the north-east corner.
zoomlevel (int): Zoom level of the tile (>=1).
filename (str): Filename for exporting the GeoDataFrame as a Shapefile.
The output starts with the North-West corner and moves first to the East, to then
jump to the row to the South, and then goes again from West to East.
Returns:
out_df (DataFrame): DataFrame with a column for the tile IDs and another for their
geometries in Well-Known Text format.
"""
print(' Gathering cell IDs...')
all_cells_ids= get_cells_in_bbox(sw_lon, sw_lat, ne_lon, ne_lat)
geom_list= []
ids_list= []
print(' Gathering geometry of cells...')
for i, id_i in enumerate(all_cells_ids):
if len(all_cells_ids)>50000:
if i%50000==0:
print(' Working on cell '+str(i+1)+' of '+str(len(all_cells_ids)))
lon_w, lat_s, lon_e, lat_n= get_coordinates_from_cell_id(id_i)
geom_list.append((shapely.geometry.Polygon([(lon_w,lat_s), (lon_e,lat_s), (lon_e,lat_n), (lon_w,lat_n)])).wkt)
ids_list.append(str(id_i))
if sw_lon > ne_lon:
print("ERROR in get_cells_in_bbox: sw_lon should be smaller than ne_lon")
return []
if sw_lat > ne_lat:
print("ERROR in get_cells_in_bbox: sw_lat should be smaller than ne_lat")
return []
if zoomlevel < 1:
print("ERROR in get_cells_in_bbox: zoomlevel needs to be equal to or larger than 1")
return []
all_tiles = mercantile.tiles(adjust_lon(sw_lon),
adjust_lat(sw_lat),
adjust_lon(ne_lon),
adjust_lat(ne_lat),
zoomlevel)
cells_ids = []
geom_list = []
print(' Gathering IDs and geometries of tiles...')
for tile in all_tiles:
cells_ids.append(mercantile.quadkey(tile))
tile_bounds = mercantile.bounds(tile)
geom_list.append((shapely.geometry.Polygon([(tile_bounds.west,
tile_bounds.south),
(tile_bounds.east,
tile_bounds.south),
(tile_bounds.east,
tile_bounds.north ),
(tile_bounds.west,
tile_bounds.north )])).wkt)
print(' Creating GeoPandas DataFrame...')
data= {'cell_id': ids_list,
data= {'cell_id': cells_ids,
'geometry': geom_list}
out_df= pd.DataFrame(data, columns = ['cell_id','geometry'])
return out_df
def get_area_polygon_on_Earth(polygon_in,units='m2'):
"""
polygon_in= Shapely polygon, defined in epsg:4326 ("normal" latitude and longitude)
units= the units of the output (default m2)
"""This function returns the area in m2 or km2 of a polygon on Earth whose geometry is
defined in epsg:4326 ("normal" latitude and longitude).
This function was written based on the code in:
https://gist.github.com/robinkraft/c6de2f988c9d3f01af3c
(general process)
and
https://stackoverflow.com/questions/51554602/how-do-i-get-the-area-of-a-geojson-polygon-with-python
https://stackoverflow.com/questions/51554602/
how-do-i-get-the-area-of-a-geojson-polygon-with-python
(improved projection, instead of epsg:3857 as suggested above)
and
https://stackoverflow.com/questions/4681737/how-to-calculate-the-area-of-a-polygon-on-the-earths-surface-using-python
(adding parameters apart from lat1 and lat2 to partial(... pyproj.Proj), which are used to define the
Albers Equal Area projection)
and
https://stackoverflow.com/questions/4681737/
how-to-calculate-the-area-of-a-polygon-on-the-earths-surface-using-python
(adding parameters apart from lat1 and lat2 to partial(... pyproj.Proj), which are used to
define the Albers Equal Area projection)
Args:
polygon_in (Shapely polygon): Polygon defined in epsg:4326.
units (str): The units of the output, can be m2 or km2 (default: m2).
Returns:
projected_area: Area of the polygon in m2 or km2.
"""
bbox= polygon_in.bounds # sw_lon, sw_lat, ne_lon, ne_lat
# Create projection from epsg:4326' to Albers Equal Area (https://proj.org/operations/projections/aea.html)
# For this it is necessary to provide as input the area of the world of interest, which is input herein from
# the bounding box of the polygon whose area we want to calculate.
proj = partial(pyproj.transform,
bbox = polygon_in.bounds # sw_lon, sw_lat, ne_lon, ne_lat
# Create projection from epsg:4326' to Albers Equal Area
# (https://proj.org/operations/projections/aea.html)
# For this it is necessary to provide as input the area of the world of interest, which is
# input herein from the bounding box of the polygon whose area we want to calculate.
inputstr = '+proj=aea +lat_1={} +lat_2={} +lat_0={} +lon_0={}'.format(bbox[1],
bbox[3],
(bbox[1]+bbox[3])/2.0,
(bbox[0]+bbox[2])/2.0)
proj = partial(pyproj.transform,
pyproj.Proj(init='epsg:4326'),
pyproj.Proj(proj='aea',lat1=bbox[1], lat2=bbox[3], lat0=(bbox[1]+bbox[3])/2.0, lon0=(bbox[0]+bbox[2])/2.0))
pyproj.Proj(inputstr))
projected_area = shapely.ops.transform(proj, polygon_in).area # in m2
if units=='km2':
projected_area= projected_area / 1000000.0
elif units!='m2':
if units == 'km2':
projected_area = projected_area / 1000000.0
elif units != 'm2':
print(' UNITS NOT RECOGNISED IN get_area_km2_polygon_on_Earth FUNCTION. ERROR!!!')
projected_area= -999.9
return projected_area
\ No newline at end of file
return projected_area
Supports Markdown
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