Source code for climada_petals.hazard.emulator.geo

This file is part of CLIMADA.

Copyright (C) 2017 ETH Zurich, CLIMADA contributors listed in AUTHORS.

CLIMADA is free software: you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free
Software Foundation, version 3.

CLIMADA 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 General Public License for more details.

You should have received a copy of the GNU General Public License along
with CLIMADA. If not, see <>.


Georegion objects for the hazard event emulator.

import logging
import numpy as np
import geopandas as gpd
import shapely.ops
import shapely.vectorized
from shapely.geometry import Polygon

from climada.hazard import Centroids
import climada.util.coordinates as u_coord
import climada_petals.hazard.emulator.const as const

LOGGER = logging.getLogger(__name__)

[docs] class HazRegion(): """Hazard region for given geo information"""
[docs] def __init__(self, extent=None, geometry=None, country=None, season=(1, 12)): """Initialize HazRegion If several arguments are passed, the spatial intersection is taken. Parameters ---------- extent : tuple (lon_min, lon_max, lat_min, lat_max), optional geometry : GeoPandas DataFrame, optional country : str or list of str, optional Countries are represented by their ISO 3166-1 alpha-3 identifiers. The keyword "all" chooses all countries (i.e., global land areas). season : pair of int, optional First and last month of hazard-specific season within this region """ self._set_geometry(extent=extent, geometry=geometry, country=country) self.geometry['const'] = 0 self.shape = self.geometry.dissolve(by='const').geometry[0] self.season = season
def _set_geometry(self, extent=None, geometry=None, country=None): self.meta = {} if extent is not None: self.meta['extent'] = extent else: extent = (-180, 180, -90, 90) lon_min, lon_max, lat_min, lat_max = extent extent_poly = gpd.GeoSeries(Polygon([ (lon_min, lat_min), (lon_min, lat_max), (lon_max, lat_max), (lon_max, lat_min) ]), crs=u_coord.NE_CRS) self.geometry = gpd.GeoDataFrame({'geometry': extent_poly}, crs=u_coord.NE_CRS) if country is not None: self.meta['country'] = country if country == "all": country = None elif not isinstance(country, list): country = [country] country_geom = u_coord.get_country_geometries(country_names=country) self.geometry = gpd.overlay(self.geometry, country_geom, how="intersection") if geometry is not None: self.meta['geometry'] = repr(geometry) self.geometry = gpd.overlay(self.geometry, geometry, how="intersection")
[docs] def centroids(self, latlon=None, res_as=360): """Return centroids in this region Parameters ---------- latlon : pair (lat, lon), optional Latitude and longitude of centroids. If not given, values are taken from CLIMADA's base grid (see `res_as`). res_as : int, optional One of 150 or 360. When `latlon` is not given, choose coordinates from centroids according to CLIMADA's base grid of given resolution in arc-seconds. Default: 360. Returns ------- centroids : climada.hazard.Centroids object """ if latlon is None: centroids = Centroids.from_base_grid(res_as=res_as) centroids.set_meta_to_lat_lon() else: centroids = Centroids.from_lat_lon(*latlon) msk = shapely.vectorized.contains(self.shape, centroids.lon, centroids = = np.arange(centroids.lon.shape[0]) return centroids
[docs] class TCRegion(HazRegion): """Hazard region with support for TC ocean basins"""
[docs] def __init__(self, tc_basin=None, season=None, **kwargs): """Initialize TCRegion The given geo information must be such that everything is contained in a single TC ocean basin. Parameters ---------- tc_basin : str TC (sub-)basin abbreviated name, such as "SIW". If not given, automatically determined from geometry and basin bounds. **kwargs : see HazRegion.__init__ """ self._set_geometry(**kwargs) self.tc_basin = None if tc_basin is not None: tc_basin_geom = get_tc_basin_geometry(tc_basin) self.geometry = gpd.overlay(self.geometry, tc_basin_geom, how="intersection") self.meta['tc_basin'] = tc_basin self.tc_basin = tc_basin self.geometry['const'] = 0 self.shape = self.geometry.dissolve(by='const').geometry[0] if self.tc_basin is None: self._determine_tc_basin() self.hemisphere = 'S' if const.TC_BASIN_GEOM_SIMPL[self.tc_basin][0][3] <= 0 else 'N' if season is None: season = const.TC_BASIN_SEASONS[self.tc_basin[:2]] self.season = season
def _determine_tc_basin(self): for basin in const.TC_SUBBASINS: basin_geom = get_tc_basin_geometry(basin) if all(basin_geom.contains(self.shape)): self.tc_basin = basin break if self.tc_basin is None: raise ValueError("Region is not contained in a single basin!") for tc_basin in const.TC_SUBBASINS[self.tc_basin]: tc_basin_geom = get_tc_basin_geometry(tc_basin) if all(tc_basin_geom.contains(self.shape)): self.tc_basin = tc_basin break"Automatically determined TC basin: %s", self.tc_basin)
[docs] def get_tc_basin_geometry(tc_basin): """Get TC (sub-)basin geometry Parameters ---------- tc_basin : str TC (sub-)basin abbreviated name, such as "SIW" or "NA". Returns ------- df : GeoPandas DataFrame """ polygons = [] for rect in const.TC_BASIN_GEOM[tc_basin]: lonmin, lonmax, latmin, latmax = rect polygons.append(Polygon([ (lonmin, latmin), (lonmin, latmax), (lonmax, latmax), (lonmax, latmin) ])) polygons = shapely.ops.unary_union(polygons) polygons = gpd.GeoSeries(polygons, crs=u_coord.NE_CRS) return gpd.GeoDataFrame({'geometry': polygons}, crs=u_coord.NE_CRS)