Commit 9b2bf257 authored by Cecilia Nievas's avatar Cecilia Nievas
Browse files

Merge branch 'feature/quadtiles' into 'master'

Adaptation of the GDE prototype code to work with quadtiles instead of 10-arcsec cells

See merge request !18
parents c8a2a0b8 79f71100
......@@ -914,7 +914,7 @@ def fill_in_admin_summary_dictionary_with_OBM(adm_summ_dict, obm_dict, obm_list_
adm_lev= obm_country_adminIDs_levels[j]
if adm_lev>-1:
country_name= gdet_serah.get_country_name_from_iso2_and_metadata_HDF5(adm_id.split('_')[0], sera_meta_path)
adm_summ_dict[adm_id]['geometry'], _= get_admin_unit_geometry_from_ID(country_name, adm_lev, adm_id.split('_')[1], countries_shp_path)
adm_summ_dict[adm_id]['geometry'], _= get_admin_unit_geometry_from_ID(country_name, adm_lev, '_'.join(adm_id.split('_')[1:]), countries_shp_path)
else:
adm_summ_dict[adm_id]['geometry']= 'NotAvailable'
return adm_summ_dict
......@@ -979,7 +979,7 @@ def fill_in_admin_summary_dictionary_with_SERA_LeftOver_or_Total_bdgs_dwells_ppl
adm_lev= list_country_adminIDs_levels[j]
if adm_lev>-1:
country_name= gdet_serah.get_country_name_from_iso2_and_metadata_HDF5(country_admin_ID.split('_')[0], sera_meta_path)
adm_summ_dict[country_admin_ID]['geometry'], _= get_admin_unit_geometry_from_ID(country_name, adm_lev, country_admin_ID.split('_')[1], countries_shp_path)
adm_summ_dict[country_admin_ID]['geometry'], _= get_admin_unit_geometry_from_ID(country_name, adm_lev, '_'.join(country_admin_ID.split('_')[1:]), countries_shp_path)
else:
adm_summ_dict[country_admin_ID]['geometry']= 'NotAvailable'
return adm_summ_dict
......@@ -1261,7 +1261,7 @@ def define_grid_cell_list_to_process(configuration_dict, script_name):
"""
out_grid_cells_def_method= configuration_dict[script_name]['method']
if out_grid_cells_def_method=='arbitrary':
out_grid_cells_list= [int(i) for i in configuration_dict[script_name]['arbitr_list'].split(', ')]
out_grid_cells_list= [i for i in configuration_dict[script_name]['arbitr_list'].split(', ')]
elif out_grid_cells_def_method=='country_in_DB':
out_grid_cells_list= gdet_psql.retrieve_cell_IDs_by_country(configuration_dict['Tiles Database']['db_tiles_name'],
configuration_dict['Tiles Database']['db_tiles_username'],
......
......@@ -27,6 +27,7 @@ Tools used by the GDE code to access/query/write to the OBM extract and cells da
"""
import psycopg2 # Documentation: http://initd.org/psycopg/docs/
import getpass
import numpy as np
import pandas as pd
from copy import deepcopy
......@@ -74,7 +75,7 @@ def retrieve_parameters_of_buildings_in_cell_raw(db_name, db_username, sch_name,
conn.set_session(autocommit=True)
cur = conn.cursor() # Open a cursor to perform database operations
aux_query_str= ','.join(parameter_names)
cur.execute('SELECT %s FROM %s.%s WHERE cell_id=%s;' %(aux_query_str, sch_name, table_name, in_cell_id))
cur.execute("SELECT %s FROM %s.%s WHERE cell_id='%s';" %(aux_query_str, sch_name, table_name, in_cell_id))
exec_result = cur.fetchall()
conn.close()
return exec_result
......@@ -617,6 +618,205 @@ def get_obm_buildings_by_broad_occupancy(db_name, db_username, sch_name, table_n
return out_broad_occus, out_num_bdgs, gem_occus, num_entries
def get_ghs_built_up_area_from_obm_tiles(grid_cell_id, db_hostname, db_name,
db_username, db_table, source_id=1):
"""This function retrieves the GHS built-up area from the OBM tiles database for a specified
zoom level 18 quadkey.
Note that a password is required to access the database. If not provided as input, this
function will prompt the user to indicate it. If the connection to the database and a cursor
are opened before calling this function, consider using the function
get_ghs_built_up_area_from_obm_tiles_opened_connection() instead.
Args:
grid_cell_id (str): Quadkey of the zoom level 18 target tile.
db_hostname (str): Name of the server where the database is located.
db_name (str): Name of the database.
db_username (str): Username for access to the database.
db_table (str): Name of the table where this information is stored.
source_id (int): Integer indicating the source ID of the GHS built-up area. Deafault: 1.
Returns:
ghs_built_up (float): Built-up area (m2) in the target tile according to the Global
Human Settlement GHSL-BUILT multitemporal dataset available in
30-m resolution at https://ghsl.jrc.ec.europa.eu/download.php.
"""
# Connect to the database:
connection = psycopg2.connect('host=%s dbname=%s user=%s password=%s'
% (db_hostname, db_name, db_username, getpass.getpass()))
connection.set_session(autocommit=True)
# Open a cursor to perform database operations
cursor = connection.cursor()
# Query the database:
sql_query = "SELECT built_area_size FROM %s WHERE (quadkey='%s' AND source_id=%s);" \
% (db_table, grid_cell_id, source_id)
cursor.execute(sql_query)
exec_result = cursor.fetchall()
connection.close()
if len(exec_result) == 0: # If the quadkey is not found the built-up area is zero
ghs_built_up = 0.0
elif len(exec_result) == 1:
ghs_built_up = exec_result[0][0]
else: # More than one entries found, this is an error
print("ERROR IN get_ghs_built_up_area_from_obm_tiles: "
"MORE THAN ONE ENTRY FOUND FOR CELL ID %s"
% (grid_cell_id))
ghs_built_up = -999.9
return ghs_built_up
def get_ghs_built_up_area_from_obm_tiles_opened_connection(grid_cell_id, open_cursor,
db_table, source_id=1):
"""This function retrieves the GHS built-up area from the OBM tiles database for a specified
zoom level 18 quadkey.
Note that this function requires that a connection to the database has been established and
a cursor opened. The function get_ghs_built_up_area_from_obm_tiles() can be used if this is
not the case.
Args:
grid_cell_id (str): Quadkey of the zoom level 18 target tile.
open_cursor: Cursor connecting to the OBM tiles database (already open).
db_table (str): Name of the table where this information is stored.
source_id (int): Integer indicating the source ID of the GHS built-up area. Deafault: 1.
Returns:
ghs_built_up (float): Built-up area (m2) in the target tile according to the Global
Human Settlement GHSL-BUILT multitemporal dataset available in
30-m resolution at https://ghsl.jrc.ec.europa.eu/download.php.
"""
# Query the database:
sql_query = "SELECT built_area_size FROM %s WHERE (quadkey='%s' AND source_id=%s);" \
% (db_table, grid_cell_id, source_id)
open_cursor.execute(sql_query)
exec_result = open_cursor.fetchall()
if len(exec_result) == 0: # If the quadkey is not found the built-up area is zero
ghs_built_up = 0.0
elif len(exec_result) == 1:
ghs_built_up = exec_result[0][0]
else: # More than one entries found, this is an error
print("ERROR IN get_ghs_built_up_area_from_obm_tiles: "
"MORE THAN ONE ENTRY FOUND FOR CELL ID %s"
% (grid_cell_id))
ghs_built_up = -999.9
return ghs_built_up
def get_completeness_from_obm_tiles(grid_cell_id, db_hostname, db_name,
db_username, db_table):
"""This function retrieves the completeness value from the OBM tiles database for a
specified zoom level 18 quadkey.
Note that a password is required to access the database. If not provided as input, this
function will prompt the user to indicate it.
Args:
grid_cell_id (str): Quadkey of the zoom level 18 target tile.
db_hostname (str): Name of the server where the database is located.
db_name (str): Name of the database.
db_username (str): Username for access to the database.
db_table (str): Name of the table where this information is stored.
Returns:
completeness (int): Completeness value in the target tile. The possible values are:
0: unknown
1: complete
2: almost complete
3: incomplete
4: undecidable
5: water
6: empty
"""
# Connect to the database:
connection = psycopg2.connect('host=%s dbname=%s user=%s password=%s'
% (db_hostname, db_name, db_username, getpass.getpass()))
connection.set_session(autocommit=True)
# Open a cursor to perform database operations
cursor = connection.cursor()
# Query the database:
sql_query = "SELECT completeness FROM %s WHERE quadkey='%s';" \
% (db_table, grid_cell_id)
cursor.execute(sql_query)
exec_result = cursor.fetchall()
connection.close()
if len(exec_result) == 0: # If the quadkey is not found, assume completeness unknown
completeness = 0
elif len(exec_result) == 1:
completeness = exec_result[0][0]
else: # More than one entries found, this is an error
print("ERROR IN get_completeness_from_obm_tiles: "
"MORE THAN ONE ENTRY FOUND FOR CELL ID %s"
% (grid_cell_id))
completeness = -999.9
return completeness
def get_completeness_from_obm_tiles_for_several(grid_cell_ids, db_hostname, db_name,
db_username, db_table):
"""This function retrieves the completeness value from the OBM tiles database for a
specified zoom level 18 quadkey.
Note that a password is required to access the database. If not provided as input, this
function will prompt the user to indicate it.
Args:
grid_cell_ids (list of str): List of quadkeys of the zoom level 18 target tiles.
db_hostname (str): Name of the server where the database is located.
db_name (str): Name of the database.
db_username (str): Username for access to the database.
db_table (str): Name of the table where this information is stored.
Returns:
completeness (list of int): List of completeness values of the target tiles. The
possible values are:
0: unknown
1: complete
2: almost complete
3: incomplete
4: undecidable
5: water
6: empty
"""
# Connect to the database:
connection = psycopg2.connect('host=%s dbname=%s user=%s password=%s'
% (db_hostname, db_name, db_username, getpass.getpass()))
connection.set_session(autocommit=True)
# Open a cursor to perform database operations
cursor = connection.cursor()
completeness = []
for grid_cell_id in grid_cell_ids:
# Query the database:
sql_query = "SELECT completeness FROM %s WHERE quadkey='%s';" \
% (db_table, grid_cell_id)
cursor.execute(sql_query)
exec_result = cursor.fetchall()
if len(exec_result) == 0: # If the quadkey is not found, assume completeness unknown
completeness.append(0)
elif len(exec_result) == 1:
completeness.append(exec_result[0][0])
else: # More than one entries found, this is an error
print("ERROR IN get_completeness_from_obm_tiles: "
"MORE THAN ONE ENTRY FOUND FOR CELL ID %s"
% (grid_cell_id))
completeness.append(-999.9)
connection.close()
return completeness
......@@ -775,9 +775,21 @@ def transform_dict_into_dataframe(dict_of_subclasses, bdg_classes):
collect_values[b,col]= dict_of_subclasses[bdgclass]['Proportions'][which[0]]
elif len(which)>1:
print('ERROR!! SOMETHING WENT WRONG IN transform_dict_into_dataframe, length='+str(len(which)))
if np.any(np.logical_or(collect_values.sum(axis=1)>1.0001,collect_values.sum(axis=1)<0.9999)):
print('WARNING IN transform_dict_into_dataframe: THE SUM PER ROW IS NOT 1.0 IN AT LEAST ONE ROW!!!')
errors_found= True
sum_per_row = collect_values.sum(axis=1)
for row in range(0,collect_values.shape[0]):
# The row does not add up to 1 but it also does not add up to zero (measured as < 1E-5):
if np.logical_or(sum_per_row[row] > 1.0001, sum_per_row[row] < 0.9999) and (sum_per_row[row] > 1E-5):
print('WARNING IN transform_dict_into_dataframe: THE SUM PER ROW IS NOT 1.0 IN AT LEAST ONE ROW!!!')
errors_found= True
elif sum_per_row[row] < 1E-5:
# If the row sums up to zero (because there are zero buildings in these admin units), replace the
# zeroes with a uniform distribution:
aux_ids = dict_of_subclasses[bdg_classes[row]]['Country_Admin_IDs'] # admin IDs for which this building class exists
uniform_proportion = 1.0/float(len(aux_ids))
for admin_id in aux_ids:
which_col = np.where(unique_involved_ids == admin_id)[0][0]
collect_values[row, which_col] = uniform_proportion
return bdg_classes, unique_involved_ids, collect_values, errors_found
......
......@@ -149,6 +149,10 @@ def check_parameters(config, section_name):
raise IOError(
"ERROR!! version_of_SERA_metadata PARAMETER MISSING FROM CONFIG FILE!!"
)
if not config.has_option("GDE_gather_SERA_and_OBM", "read_completeness"):
raise IOError(
"ERROR!! read_completeness PARAMETER MISSING FROM CONFIG FILE!!"
)
if not config.has_option("GDE_gather_SERA_and_OBM", "print_screen_during_run"):
raise IOError(
"ERROR!! print_screen_during_run PARAMETER MISSING FROM CONFIG FILE!!"
......@@ -291,11 +295,36 @@ def check_parameters(config, section_name):
elif section_name == "SERA_mapping_admin_units_to_cells":
if not config.has_option("SERA_mapping_admin_units_to_cells", "countries"):
raise IOError("ERROR!! countries PARAMETER MISSING FROM CONFIG FILE!!")
elif section_name == "SERA_mapping_admin_units_to_cells_add_GHS":
elif section_name == "SERA_mapping_admin_units_to_cells_add_GHS_from_HDF5":
if not config.has_option(
"SERA_mapping_admin_units_to_cells_add_GHS", "ghs_input_filename"
"SERA_mapping_admin_units_to_cells_add_GHS_from_HDF5", "ghs_input_filename"
):
raise IOError("ERROR!! ghs_input_filename PARAMETER MISSING FROM CONFIG FILE!!")
elif section_name == "OBM tiles":
if not config.has_option(
"OBM tiles", "obm_tiles_db_hostname"
):
raise IOError("ERROR!! obm_tiles_db_hostname PARAMETER MISSING FROM CONFIG FILE!!")
if not config.has_option(
"OBM tiles", "obm_tiles_db_name"
):
raise IOError("ERROR!! obm_tiles_db_name PARAMETER MISSING FROM CONFIG FILE!!")
if not config.has_option(
"OBM tiles", "obm_tiles_db_username"
):
raise IOError("ERROR!! obm_tiles_db_username PARAMETER MISSING FROM CONFIG FILE!!")
if not config.has_option(
"OBM tiles", "obm_tiles_db_table_built_up"
):
raise IOError("ERROR!! obm_tiles_db_table_built_up PARAMETER MISSING FROM CONFIG FILE!!")
if not config.has_option(
"OBM tiles", "obm_tiles_db_table_completeness"
):
raise IOError("ERROR!! obm_tiles_db_table_completeness PARAMETER MISSING FROM CONFIG FILE!!")
if not config.has_option(
"OBM tiles", "obm_tiles_source_built_up_id"
):
raise IOError("ERROR!! obm_tiles_source_id PARAMETER MISSING FROM CONFIG FILE!!")
elif section_name == "SERA_mapping_admin_units_to_cells_add_Sat":
if not config.has_option("SERA_mapping_admin_units_to_cells_add_Sat", "apply_model"):
raise IOError("ERROR!! apply_model PARAMETER MISSING FROM CONFIG FILE!!")
......
"""
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):