GDE_TOOLS_world_grid.py 13.9 KB
Newer Older
1
"""
2
3
Copyright (C) 2021
  Helmholtz-Zentrum Potsdam Deutsches GeoForschungsZentrum GFZ
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

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/>.


Global Dynamic Exposure Model
Helmholtz Centre Potsdam
GFZ German Research Centre for Geosciences
Section 2.6: Seismic Hazard and Risk Dynamics

GDE_TOOLS_world_grid
====================
26
27
28
29
30
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.

31
32
33
34
35
36
37
38
"""

import numpy as np
import geopandas as gpd
import pandas as pd
import shapely
import pyproj
from functools import partial
39
40
import mercantile

41
42

def adjust_lon(coord_lon):
43
44
45
46
47
48
49
    """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].    
50
    """
51
52
53

    if coord_lon<0.0:
        return max(coord_lon, -180.0)
54
    else:
55
56
        return min(coord_lon, 180.0)

57
58
    
def adjust_lat(coord_lat):
59
60
61
62
63
64
65
    """ 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].       
66
    """
67
68
69

    if coord_lat<0.0:
        return max(coord_lat, -85.0511)
70
    else:
71
        return min(coord_lat, 85.0511)
72

73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89

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.
90
    """
91
92
93
94
95
96
97
98
99

    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))

100
101
    return cell_id

102
103
104
105
106
107
108
109
110
111
112
113
114
115

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)
116
    """
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138

    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))                  

139
140
    return cells_ids

141

142
def get_coordinates_from_cell_id(cell_id):
143
144
145
146
147
148
149
150
151
152
153
    """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.
154
    """
155
156
157
158
159
160
161

    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       

162
163
164
165
    return lon_w, lat_s, lon_e, lat_n


def get_cells_geometry_as_geopandas(cell_id_list, filename=None):
166
167
168
169
170
171
172
173
174
    """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.
175
    """
176
177

    geom_list = []
178
179
    for cell_id in cell_id_list:
        lon_w, lat_s, lon_e, lat_n= get_coordinates_from_cell_id(cell_id)
180
181
182
183
        geom_list.append(shapely.geometry.Polygon([(lon_w, lat_s),
                                                   (lon_e, lat_s),
                                                   (lon_e, lat_n),
                                                   (lon_w, lat_n)]))
184
185
186
187
188
189
190
191
    geometry = gpd.GeoSeries(geom_list) 
    data= {'cell_id': cell_id_list,
           'geometry': geometry}
    out_gdf= gpd.GeoDataFrame(data, columns = ['cell_id','geometry'])
    out_gdf.crs = {'init' :'epsg:4326'}
    if filename:
        print('      Saving to GeoJSON...')
        out_gdf.to_file(filename,driver='GeoJSON')
192

193
194
195
196
    return out_gdf    


def get_geometry_WKT_of_cell(cell_id):
197
198
199
    """This function returns the Well-Known Text of the geometry of the tile with cell_id.
    """

200
    lon_w, lat_s, lon_e, lat_n= get_coordinates_from_cell_id(cell_id)
201
202
203
204
205
    out_wkt= (shapely.geometry.Polygon([(lon_w, lat_s),
                                        (lon_e, lat_s),
                                        (lon_e, lat_n),
                                        (lon_w, lat_n)])).wkt

206
    return out_wkt
207
208


209
def get_cells_geometry_as_pandas_with_WKT(cell_id_list):
210
211
212
213
214
215
216
217
218
    """ 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.
219
    """
220
221

    geom_list = []
222
223
    for cell_id in cell_id_list:
        lon_w, lat_s, lon_e, lat_n= get_coordinates_from_cell_id(cell_id)
224
225
226
227
228
        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,
229
           'geometry': geom_list}
230
231
    out_df = pd.DataFrame(data, columns = ['cell_id','geometry'])

232
233
    return out_df  

234
235
236
237
238
239
240
241
242
243
244
245
246

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.
247
    
248
249
250
    Returns:
        out_gdf (GeoDataFrame): GeoDataFrame with the tile IDs and another with their
                                geometries.
251
    """
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280

    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 )]))

281
282
    print('      Creating GeoPandas DataFrame...')
    geometry = gpd.GeoSeries(geom_list) 
283
    data= {'cell_id': cells_ids,
284
285
286
           'geometry': geometry}
    out_gdf= gpd.GeoDataFrame(data, columns = ['cell_id','geometry'])
    out_gdf.crs = {'init' :'epsg:4326'}
287

288
289
290
    if filename:
        print('      Saving to ShapeFile...')
        out_gdf.to_file(filename,driver='ESRI Shapefile')
291

292
293
294
    return out_gdf    

    
295
296
297
298
299
300
301
302
303
304
305
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.
306
    
307
308
309
    Returns:
        out_df (DataFrame): DataFrame with a column for the tile IDs and another for their
                            geometries in Well-Known Text format.
310
    """
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343

    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)

344
    print('      Creating GeoPandas DataFrame...')
345
    data= {'cell_id': cells_ids,
346
347
           'geometry': geom_list}
    out_df= pd.DataFrame(data, columns = ['cell_id','geometry'])
348

349
350
351
    return out_df    
    

352
353
354
def get_area_polygon_on_Earth(input_polygon, units='m2'):
    """This function returns the area (in m2 or km2) of a polygon on Earth calculated using the
    Albers equal area projection. The polygon's geometry is assumed to be defined in ESPG:4326.
355
356

    Args:
357
        input_polygon (Shapely polygon): Polygon defined in epsg:4326.
358
359
360
361
        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.
362
    """
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
    
    # 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)

    projected_area = geometry.area  # in m2
379
380
381
382

    if units == 'km2':
        projected_area = projected_area / 1000000.0
    elif units != 'm2':
383
384
        print('      UNITS NOT RECOGNISED IN get_area_km2_polygon_on_Earth FUNCTION. ERROR!!!')
        projected_area= -999.9
385

386
    return projected_area