"""
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 <https://www.gnu.org/licenses/>.
---
Define Centroids class.
"""
import copy
from deprecation import deprecated
import logging
from pathlib import Path
from typing import Any, Literal, Union
import warnings
import h5py
import cartopy.crs as ccrs
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pyproj.crs.crs import CRS
import rasterio
from shapely.geometry.point import Point
from climada.util.constants import DEF_CRS
import climada.util.coordinates as u_coord
import climada.util.plot as u_plot
__all__ = ['Centroids']
PROJ_CEA = CRS.from_user_input({'proj': 'cea'})
LOGGER = logging.getLogger(__name__)
DEF_SHEET_NAME = 'centroids'
EXP_SPECIFIC_COLS = [
'value',
'impf_',
'centr_',
'cover',
'deductible'
]
[docs]
class Centroids():
"""Contains vector centroids as a GeoDataFrame
Attributes
----------
lat : np.array
Latitudinal coordinates in the specified CRS (can be any unit).
lon : np.array
Longitudinal coordinates in the specified CRS (can be any unit).
crs : pyproj.CRS
Coordinate reference system. Default: EPSG:4326 (WGS84)
region_id : np.array, optional
Numeric country (or region) codes. Default: None
on_land : np.array, optional
Boolean array indicating on land (True) or off shore (False). Default: None
"""
[docs]
def __init__(
self,
*,
lat: Union[np.ndarray, list[float]],
lon: Union[np.ndarray, list[float]],
crs: Any = DEF_CRS,
region_id: Union[Literal["country"], None, np.ndarray, list[float]] = None,
on_land: Union[Literal["natural_earth"], None, np.ndarray, list[bool]] = None,
**kwargs,
):
"""Initialization
Parameters
----------
lat : np.array
Latitudinal coordinates in the specified CRS (can be any unit).
lon : np.array
Longitudinal coordinates in the specified CRS (can be any unit).
crs : str or anything accepted by pyproj.CRS.from_user_input()
Coordinate reference system. Default: EPSG:4326 (WGS84)
region_id : np.array or str, optional
Array of numeric country (or region) codes. If the special value "country" is given
admin-0 codes are automatically assigned. Default: None
on_land : np.array or str, optional
Boolean array indicating on land (True) or off shore (False). If the special value
"natural_earth" is given, the property is automatically determined from NaturalEarth
shapes. Default: None
kwargs : dict
Additional columns with data to store in the internal GeoDataFrame (gdf attribute).
"""
self.gdf = gpd.GeoDataFrame(
data={
'geometry': gpd.points_from_xy(lon, lat, crs=crs),
'region_id': region_id,
'on_land': on_land,
**kwargs,
}
)
if isinstance(region_id, str):
LOGGER.info('Setting region id to %s level.', region_id)
self.set_region_id(level=region_id, overwrite=True)
if isinstance(on_land, str):
LOGGER.info('Setting on land from %s source.', on_land)
self.set_on_land(source=on_land, overwrite=True)
@property
def lat(self):
""" Return latitudes """
return self.gdf.geometry.y.values
@property
def lon(self):
""" Return longitudes """
return self.gdf.geometry.x.values
@property
def geometry(self):
""" Return the geometry """
return self.gdf['geometry']
@property
def on_land(self):
""" Get the on_land property """
if self.gdf.on_land.isna().all():
return None
return self.gdf['on_land'].values
@property
def region_id(self):
""" Get the assigned region_id """
if self.gdf.region_id.isna().all():
return None
return self.gdf['region_id'].values
@property
def crs(self):
""" Get the crs"""
return self.gdf.crs
@property
def size(self):
"""Get size (number of lat/lon pairs)"""
return self.gdf.shape[0]
@property
def shape(self):
"""Get shape [lat, lon] assuming rastered data."""
return (np.unique(self.lat).size, np.unique(self.lon).size)
@property
def total_bounds(self):
"""Get total bounds (minx, miny, maxx, maxy)."""
return self.gdf.total_bounds
@property
def coord(self):
"""Get [lat, lon] array."""
return np.stack([self.lat, self.lon], axis=1)
def __eq__(self, other):
""" dunder method for Centroids comparison.
returns True if two centroids equal, False otherwise
Parameters
----------
other : Centroids
object to compare with
Returns
-------
bool
"""
if not u_coord.equal_crs(self.crs, other.crs):
return False
try:
pd.testing.assert_frame_equal(self.gdf, other.gdf, check_like=True)
return True
except AssertionError:
return False
[docs]
def to_default_crs(self, inplace=True):
"""Project the current centroids to the default CRS (epsg4326)
Parameters
----------
inplace: bool
if True, modifies the centroids in place.
if False, return projected centroids object.
Default is True.
Returns
-------
Centroids or None (if inplace is True)
"""
return self.to_crs(DEF_CRS, inplace=inplace)
[docs]
def to_crs(self, crs, inplace=False):
""" Project the current centroids to the desired crs
Parameters
----------
crs : str
coordinate reference system
inplace: bool, default False
if True, modifies the centroids in place.
if False, returns a copy.
Returns
-------
Centroids or None (if inplace is True)
"""
if inplace:
self.gdf.to_crs(crs, inplace=True)
return None
return Centroids.from_geodataframe(self.gdf.to_crs(crs))
[docs]
@classmethod
def from_geodataframe(cls, gdf):
"""Initialize centroids from a geodataframe
Parameters
----------
gdf : GeoDataFrame
Input geodataframe with centroids as points
in the geometry column. All other columns are
attached to the centroids geodataframe.
Returns
-------
Centroids
Centroids built from the geodataframe.
Raises
------
ValueError
"""
if (gdf.geom_type != 'Point').any():
raise ValueError(
'The inpute geodataframe contains geometries that are not points.'
)
# Don't forget to make a copy!!
# This is a bit ugly, but avoids to recompute the geometries
# in the init. For large datasets this saves computation time
centroids = cls(lat=[1], lon=[1]) #make "empty" centroids
centroids.gdf = gdf.copy(deep=True)
if gdf.crs is None:
centroids.gdf.set_crs(DEF_CRS, inplace=True)
return centroids
[docs]
@classmethod
def from_exposures(cls, exposures):
"""Generate centroids from the locations of exposures.
The properties "region_id" and "on_land" are also extracted from the Exposures object if
available. The columns "value", "impf_*", "centr_*", "cover", and "deductible" are not
used.
Parameters
----------
exposures : Exposures
Exposures from which to take the centroids locations (as well as region_id and on_land
if available).
Returns
-------
Centroids
Raises
------
ValueError
"""
col_names = [
column for column in exposures.gdf.columns
if not any(pattern in column for pattern in EXP_SPECIFIC_COLS)
]
# Legacy behaviour
# Exposures can be without geometry column
#TODO: remove once exposures is real geodataframe with geometry.
if 'geometry' in exposures.gdf.columns:
gdf = exposures.gdf[col_names]
return cls.from_geodataframe(gdf)
if 'latitude' in exposures.gdf.columns and 'longitude' in exposures.gdf.columns:
gdf = exposures.gdf[col_names]
return cls(
lat=exposures.gdf['latitude'],
lon=exposures.gdf['longitude'],
crs=exposures.crs,
**dict(gdf.items()),
)
raise ValueError(
"The given exposures object has no coordinates information."
"The exposures' GeoDataFrame must have either point geometries"
" or latitude and longitude values."
)
[docs]
@classmethod
def from_pnt_bounds(cls, points_bounds, res, crs=DEF_CRS):
"""Create Centroids object from coordinate bounds and resolution.
The result contains all points from a regular raster with the given resolution and CRS,
covering the given bounds. Note that the raster bounds are larger than the points' bounds
by res/2.
Parameters
----------
points_bounds : tuple
The bounds (lon_min, lat_min, lon_max, lat_max) of the point coordinates.
res : float
The desired resolution in same units as `points_bounds`.
crs : dict() or rasterio.crs.CRS, optional
Coordinate reference system. Default: DEF_CRS
Returns
-------
Centroids
"""
height, width, transform = u_coord.pts_to_raster_meta(points_bounds, (res, -res))
return cls.from_meta({
"crs": crs,
"width": width,
"height": height,
"transform": transform,
})
[docs]
def append(self, centr):
"""Append Centroids
Note that the result might contain duplicate points if the object to append has an overlap
with the current object.
Parameters
----------
centr : Centroids
Centroids to append. The centroids need to have the same CRS.
Raises
------
ValueError
See Also
--------
union : Union of Centroid objects.
remove_duplicate_points : Remove duplicate points in a Centroids object.
"""
if not u_coord.equal_crs(self.crs, centr.crs):
raise ValueError(
"The centroids have different Coordinate-Reference-Systems (CRS)"
)
self.gdf = pd.concat([self.gdf, centr.gdf])
[docs]
def union(self, *others):
"""Create the union of Centroids objects
All centroids must have the same CRS. Points that are contained in more than one of the
Centroids objects will only be contained once (i.e. duplicates are removed).
Parameters
----------
others : list of Centroids
Centroids contributing to the union.
Returns
-------
centroids : Centroids
Centroids object containing the union of all Centroids.
"""
centroids = copy.deepcopy(self)
for cent in others:
centroids.append(cent)
return centroids.remove_duplicate_points()
[docs]
def remove_duplicate_points(self):
"""Return a copy of centroids with duplicate points removed
Parameters
----------
centr : Centroids
Centroids with or without duplicate points
Returns
-------
centroids : Centroids
A new Centroids object that contains a subselection of the original centroids without
duplicates. Note that a copy is returned even if there were no duplicates.
"""
return self.from_geodataframe(self.gdf.drop_duplicates(subset=["geometry"]))
[docs]
def select(self, reg_id=None, extent=None, sel_cen=None):
"""Return new Centroids object containing points following certain criteria
It is currently possible to filter by region (reg_id), by geographical extent (extent), or
by an explicit list of indices/a mask (sel_cen). If more than one criterion is given, all
of them must be satisfied for a point to be included in the selection.
Parameters
----------
reg_id : int or list of int, optional
Numeric ID (or IDs) of the region (or regions) to restrict to, according to the values
in the region_id property. Default: None
extent : tuple, optional
The geographical extent (min_lon, max_lon, min_lat, max_lat) to restrict to, including
the boundary. If the value for min_lon is greater than max_lon, the extent is
interpreted to cross the antimeridian ([max_lon, 180] and [-180, min_lon]).
Default: None
sel_cen : np.ndarray of int or bool, optional
Boolean mask, or list of indices to restrict to. Default: None
Returns
-------
centroids : Centroids
Sub-selection of this object
"""
sel_cen_bool = sel_cen
if sel_cen is not None and sel_cen.dtype.kind == 'i':
# if needed, convert indices to bool
sel_cen_bool = np.zeros(self.size, dtype=bool)
sel_cen_bool[np.unique(sel_cen)] = True
sel_cen_mask = self.select_mask(sel_cen=sel_cen_bool, reg_id=reg_id, extent=extent)
return Centroids.from_geodataframe(self.gdf.iloc[sel_cen_mask])
[docs]
def select_mask(self, sel_cen=None, reg_id=None, extent=None):
"""Create mask of selected centroids
Parameters
----------
sel_cen: np.ndarray of bool, optional
Boolean array, with size matching the number of centroids. Default: None
reg_id : int or list of int, optional
Numeric ID (or IDs) of the region (or regions) to restrict to, according to the values
in the region_id property. Default: None
extent : tuple, optional
The geographical extent (min_lon, max_lon, min_lat, max_lat) to restrict to, including
the boundary. If the value for min_lon is greater than lon_max, the extent is
interpreted to cross the antimeridian ([lon_max, 180] and [-180, lon_min]).
Default: None
Returns
-------
sel_cen : np.ndarray of bool
Boolean array (mask) with value True for centroids in selection.
"""
if sel_cen is None:
sel_cen = np.ones(self.size, dtype=bool)
if reg_id is not None:
sel_cen &= np.isin(self.region_id, reg_id)
if extent is not None:
lon_min, lon_max, lat_min, lat_max = extent
lon_max += 360 if lon_min > lon_max else 0
lon_normalized = u_coord.lon_normalize(
self.lon.copy(), center=0.5 * (lon_min + lon_max))
sel_cen &= (
(lon_normalized >= lon_min) & (lon_normalized <= lon_max) &
(self.lat >= lat_min) & (self.lat <= lat_max)
)
return sel_cen
#TODO replace with nice GeoDataFrame util plot method.
[docs]
def plot(self, axis=None, figsize=(9, 13), **kwargs):
"""Plot centroids scatter points over earth
Parameters
----------
axis : matplotlib.axes._subplots.AxesSubplot, optional
axis to use
figsize: (float, float), optional
figure size for plt.subplots
The default is (9, 13)
kwargs : optional
arguments for scatter matplotlib function
Returns
-------
axis : matplotlib.axes._subplots.AxesSubplot
"""
pad = np.abs(u_coord.get_resolution(self.lat, self.lon)).min()
proj_data, _ = u_plot.get_transformation(self.crs)
proj_plot = proj_data
if isinstance(proj_data, ccrs.PlateCarree):
# use different projections for plot and data to shift the central lon in the plot
xmin, ymin, xmax, ymax = u_coord.latlon_bounds(self.lat, self.lon, buffer=pad)
proj_plot = ccrs.PlateCarree(central_longitude=0.5 * (xmin + xmax))
else:
xmin, ymin, xmax, ymax = (self.lon.min() - pad, self.lat.min() - pad,
self.lon.max() + pad, self.lat.max() + pad)
if not axis:
_fig, axis, _fontsize = u_plot.make_map(proj=proj_plot, figsize=figsize)
axis.set_extent((xmin, xmax, ymin, ymax), crs=proj_data)
u_plot.add_shapes(axis)
axis.scatter(self.lon, self.lat, transform=proj_data, **kwargs)
plt.tight_layout()
return axis
[docs]
def set_region_id(self, level='country', overwrite=False):
"""Set region_id as country ISO numeric code attribute for every pixel or point.
Parameters
----------
level: str, optional
The admin level on which to assign centroids. Currently only 'country' (admin0) is
implemented. Default: 'country'
overwrite : bool, optional
If True, overwrite the existing region_id information. If False, region_id is set
only if region_id is missing (None). Default: False
"""
if overwrite or self.region_id is None:
LOGGER.debug('Setting region_id %s points.', str(self.size))
if level == 'country':
ne_geom = self._ne_crs_geom()
self.gdf['region_id'] = u_coord.get_country_code(
ne_geom.y.values, ne_geom.x.values,
)
else:
raise NotImplementedError(
'The region id can only be assigned for countries so far'
)
[docs]
def set_on_land(self, source='natural_earth', overwrite=False):
"""Set on_land attribute for every pixel or point.
Parameters
----------
source: str, optional
The source of the on-land information. Currently, only 'natural_earth' (based on shapes
from NaturalEarth, https://www.naturalearthdata.com/) is implemented.
Default: 'natural_earth'.
overwrite : bool, optional
If True, overwrite the existing on_land information. If False, on_land is set
only if on_land is missing (None). Default: False
"""
if overwrite or self.on_land is None:
LOGGER.debug('Setting on_land %s points.', str(self.lat.size))
if source=='natural_earth':
ne_geom = self._ne_crs_geom()
self.gdf['on_land'] = u_coord.coord_on_land(
ne_geom.y.values, ne_geom.x.values
)
else:
raise NotImplementedError(
'The on land variables can only be automatically assigned using natural earth.'
)
[docs]
def get_pixel_shapes(self, res=None, **kwargs):
"""Create a GeoSeries of the quadratic pixel shapes at the centroid locations
Note that this assumes that the centroids define a regular grid of pixels.
Parameters
----------
res : float, optional
The resolution of the regular grid the pixels are taken from. If not given, it is
estimated using climada.util.coordinates.get_resolution. Default: None
kwargs : optional
Additional keyword arguments are passed to climada.util.coordinates.get_resolution.
Returns
-------
GeoSeries
See also
--------
climada.util.coordinates.get_resolution
"""
if res is None:
res = np.abs(u_coord.get_resolution(self.lat, self.lon, **kwargs)).min()
geom = self.geometry.copy()
# unset CRS to avoid warnings about geographic CRS when using `GeoSeries.buffer`
geom.crs = None
return geom.buffer(
# resolution=1, cap_style=3: squared buffers
# https://shapely.readthedocs.io/en/latest/manual.html#object.buffer
distance=res / 2, resolution=1, cap_style=3,
# reset CRS (see above)
).set_crs(self.crs)
[docs]
def get_area_pixel(self, min_resol=1.0e-8):
"""Compute the area per centroid in the CEA projection
Note that this assumes that the centroids define a regular grid of pixels (area in m²).
Parameters
----------
min_resol : float, optional
When estimating the grid resolution, use this as the minimum resolution in lat and lon.
It is passed to climada.util.coordinates.get_resolution. Default: 1.0e-8
Returns
-------
areapixels : np.array
Area of each pixel in square meters.
"""
LOGGER.debug('Computing pixel area for %d centroids.', self.size)
xy_pixels = self.get_pixel_shapes(min_resol=min_resol)
if PROJ_CEA != xy_pixels.crs:
xy_pixels = xy_pixels.to_crs(crs={'proj': 'cea'})
return xy_pixels.area.values
[docs]
def get_closest_point(self, x_lon, y_lat):
"""Returns closest centroid and its index to a given point.
Parameters
----------
x_lon : float
Longitudinal (x) coordinate.
y_lat : float
Latitudinal (y) coordinate.
Returns
-------
x_close : float
x-coordinate (longitude) of closest centroid.
y_close : float
y-coordinate (latitude) of closest centroid.
idx_close : int
Index of centroid in internal ordering of centroids.
"""
close_idx = self.geometry.distance(Point(x_lon, y_lat)).values.argmin()
return self.lon[close_idx], self.lat[close_idx], close_idx
# NOT REALLY AN ELEVATION FUNCTION, JUST READ RASTER
[docs]
def get_elevation(self, topo_path):
"""Return elevation attribute for every pixel or point in meters.
Parameters
----------
topo_path : str
Path to a raster file containing gridded elevation data.
Returns
-------
values : np.array of shape (npoints,)
Interpolated elevation values from raster file for each given coordinate point.
"""
return u_coord.read_raster_sample(topo_path, self.lat, self.lon)
[docs]
def get_dist_coast(self, signed=False, precomputed=True):
"""Get dist_coast attribute for every pixel or point in meters.
The distances are read from a raster file containing precomputed distances (from NASA) at
0.01 degree (approximately 1 km) resolution.
Parameters
----------
signed : bool, optional
If True, use signed distances (positive off shore and negative on land). Default: False
precomputed : bool, optional
Whether distances should be read from a pre-computed raster (True) or computed
on-the-fly (False). Default: True.
.. deprecated:: 5.0
Argument is ignored, because distances are not computed on-the-fly anymore.
Returns
-------
dist : np.array
(Signed) distance to coast in meters.
"""
if not precomputed:
LOGGER.warning(
"The `precomputed` argument is deprecated and will be removed in the future"
" because `get_dist_coast` always uses precomputed distances."
)
ne_geom = self._ne_crs_geom()
return u_coord.dist_to_coast_nasa(
ne_geom.y.values, ne_geom.x.values, highres=True, signed=signed,
)
##
# I/O methods
##
[docs]
@classmethod
def from_raster_file(cls, file_name, src_crs=None, window=None, geometry=None,
dst_crs=None, transform=None, width=None, height=None,
resampling=rasterio.warp.Resampling.nearest, return_meta=False):
"""Create a new Centroids object from a raster file
Select region using window or geometry. Reproject input by providing
dst_crs and/or (transform, width, height).
Parameters
----------
file_name : str
path of the file
src_crs : crs, optional
source CRS. Provide it if error without it.
window : rasterio.windows.Window, optional
window to read
geometry : list of shapely.geometry, optional
consider pixels only within these shapes
dst_crs : crs, optional
reproject to given crs
transform : rasterio.Affine
affine transformation to apply
wdith : float
number of lons for transform
height : float
number of lats for transform
resampling : rasterio.warp.Resampling optional
resampling function used for reprojection to dst_crs,
default: nearest
return_meta : bool, optional
default: False
Returns
-------
centr : Centroids
Centroids according to the given raster file
meta : dict, optional if return_meta is True
Raster meta (height, width, transform, crs).
"""
meta, _ = u_coord.read_raster(
file_name, [1], src_crs, window, geometry, dst_crs,
transform, width, height, resampling,
)
centr = cls.from_meta(meta)
return (centr, meta) if return_meta else centr
[docs]
@classmethod
def from_vector_file(cls, file_name, dst_crs=None):
"""Create Centroids object from vector file (any format supported by fiona).
Parameters
----------
file_name : str
vector file with format supported by fiona and 'geometry' field.
dst_crs : crs, optional
reproject to given crs
If no crs is given in the file, simply sets the crs.
Returns
-------
centr : Centroids
Centroids with points according to the given vector file
"""
centroids = cls.from_geodataframe(gpd.read_file(file_name))
if dst_crs is not None:
if centroids.crs:
centroids.to_crs(dst_crs, inplace=True)
else:
centroids.gdf.set_crs(dst_crs, inplace=True)
return centroids
[docs]
@classmethod
def from_csv(cls, file_path, **kwargs):
"""Generate centroids from a CSV file with column names in var_names.
Parameters
----------
file_path : str
path to CSV file to be read
kwargs : dict
Additional keyword arguments to pass on to pandas.read_csv.
Returns
-------
Centroids
"""
return cls._from_dataframe(pd.read_csv(file_path, **kwargs))
[docs]
def write_csv(self, file_path):
"""Save centroids as CSV file
Parameters
----------
file_path : str, Path
absolute or relative file path and name to write to
"""
file_path = Path(file_path).with_suffix('.csv')
LOGGER.info('Writing %s', file_path)
self._centroids_to_dataframe().to_csv(file_path, index=False)
[docs]
@classmethod
def from_excel(cls, file_path, sheet_name=None):
"""Generate a new centroids object from an excel file with column names in var_names.
Parameters
----------
file_path : str
absolute or relative file path
sheet_name : str, optional
name of sheet in excel file containing centroid information
Default: "centroids"
Returns
-------
centr : Centroids
Centroids with data from the given excel file
"""
if sheet_name is None:
sheet_name = 'centroids'
df = pd.read_excel(file_path, sheet_name)
return cls._from_dataframe(df)
[docs]
def write_excel(self, file_path):
"""Save centroids as excel file
Parameters
----------
file_path : str, Path
absolute or relative file path and name to write to
"""
file_path = Path(file_path).with_suffix('.xlsx')
LOGGER.info('Writing %s', file_path)
self._centroids_to_dataframe().to_excel(
file_path, sheet_name=DEF_SHEET_NAME, index=False,
)
[docs]
def write_hdf5(self, file_name, mode='w'):
"""Write data frame and metadata in hdf5 format
Parameters
----------
file_name : str
(path and) file name to write to.
"""
LOGGER.info('Writing %s', file_name)
store = pd.HDFStore(file_name, mode=mode)
pandas_df = pd.DataFrame(self.gdf)
for col in pandas_df.columns:
if str(pandas_df[col].dtype) == "geometry":
pandas_df[col] = np.asarray(self.gdf[col])
# Avoid pandas PerformanceWarning when writing HDF5 data
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=pd.errors.PerformanceWarning)
# Write dataframe
store.put('centroids', pandas_df)
store.get_storer('centroids').attrs.metadata = {
'crs': CRS.from_user_input(self.crs).to_wkt()
}
store.close()
[docs]
@classmethod
def from_hdf5(cls, file_name):
"""Create a centroids object from a HDF5 file.
Parameters
----------
file_data : str or h5
If string, path to read data. If h5 object, the datasets will be read from there.
Returns
-------
centr : Centroids
Centroids with data from the given file
Raises
------
FileNotFoundError
"""
if not Path(file_name).is_file():
raise FileNotFoundError(str(file_name))
try:
with pd.HDFStore(file_name, mode='r') as store:
metadata = store.get_storer('centroids').attrs.metadata
# in previous versions of CLIMADA and/or geopandas,
# the CRS was stored in '_crs'/'crs'
crs = metadata.get('crs')
gdf = gpd.GeoDataFrame(store['centroids'], crs=crs)
except TypeError:
with h5py.File(file_name, 'r') as data:
gdf = cls._gdf_from_legacy_hdf5(data.get('centroids'))
except KeyError:
with h5py.File(file_name, 'r') as data:
gdf = cls._gdf_from_legacy_hdf5(data)
return cls.from_geodataframe(gdf)
##
# Private methods
##
@classmethod
def _from_dataframe(cls, df):
if 'crs' in df.columns:
crs = df['crs'].iloc[0]
else:
LOGGER.info("No 'crs' column provided in file, setting CRS to WGS84 default.")
crs = DEF_CRS
extra_values = {
col: df[col]
for col in df.columns
if col not in ['lat', 'lon', 'crs']
}
return cls(lat=df['lat'], lon=df['lon'], **extra_values, crs=crs)
@staticmethod
def _gdf_from_legacy_hdf5(data):
crs = DEF_CRS
if data.get('crs'):
crs = u_coord.to_crs_user_input(data.get('crs')[0])
if data.get('lat') and data.get('lat').size:
latitude = np.array(data.get('lat'))
longitude = np.array(data.get('lon'))
elif data.get('latitude') and data.get('latitude').size:
latitude = np.array(data.get('latitude'))
longitude = np.array(data.get('longitude'))
else:
centr_meta = data.get('meta')
meta = dict()
meta['crs'] = crs
for key, value in centr_meta.items():
if key != 'transform':
meta[key] = value[0]
else:
meta[key] = rasterio.Affine(*value)
latitude, longitude = _meta_to_lat_lon(meta)
extra_values = {}
for centr_name in data.keys():
if centr_name not in ('crs', 'lat', 'lon', 'meta', 'latitude', 'longitude'):
values = np.array(data.get(centr_name))
if latitude.size != 0 and values.size != 0 :
extra_values[centr_name] = values
return gpd.GeoDataFrame(
extra_values,
geometry=gpd.points_from_xy(x=longitude, y=latitude, crs=crs),
)
@classmethod
def _legacy_from_excel(cls, file_name, var_names):
LOGGER.info('Reading %s', file_name)
try:
df = pd.read_excel(file_name, var_names['sheet_name'])
df = df.rename(columns=var_names['col_name'])
except KeyError as err:
raise KeyError("Not existing variable: %s" % str(err)) from err
return cls._from_dataframe(df)
def _centroids_to_dataframe(self):
"""Create dataframe from Centroids object to facilitate
saving in different file formats.
Returns
-------
df : DataFrame
"""
df = pd.DataFrame(self.gdf)
df['lon'] = self.gdf['geometry'].x
df['lat'] = self.gdf['geometry'].y
df['crs'] = CRS.from_user_input(self.crs).to_wkt()
df = df.drop(['geometry'], axis=1)
return df
def _ne_crs_geom(self):
"""Return `geometry` attribute in the CRS of Natural Earth.
Returns
-------
geo : gpd.GeoSeries
"""
if u_coord.equal_crs(self.gdf.crs, u_coord.NE_CRS):
return self.gdf.geometry
return self.to_crs(u_coord.NE_CRS, inplace=False).geometry
##
# Deprecated methods
##
[docs]
@classmethod
@deprecated(details="Reading Centroids data from matlab files is not supported anymore."
"This method has been removed with climada 5.0")
def from_mat(cls, file_name, var_names=None):
"""Reading Centroids data from matlab files is not supported anymore.
This method has been removed with climada 5.0"""
raise NotImplementedError("You are suggested to use an old version of climada (<=4.*) and"
" convert the file to hdf5 format.")
[docs]
@staticmethod
@deprecated(details="This method has been removed with climada 5.0")
def from_base_grid(land=False, res_as=360, base_file=None):
"""This method has been removed with climada 5.0"""
raise NotImplementedError("Create the Centroids from a custom base file or from Natural"
" Earth (files are available in Climada, look up ``climada.util"
".constants.NATEARTH_CENTROIDS`` for their location)")
[docs]
@classmethod
@deprecated(details="This method will be removed in a future version."
" Simply use the constructor instead.")
def from_lat_lon(cls, lat, lon, crs="EPSG:4326"):
"""deprecated, use the constructor instead"""
return cls(lat=lat, lon=lon, crs=crs)
[docs]
@deprecated(details="This method is futile and will be removed in a future version."
" `Centroids.get_area_pixel` can be run without initialization.")
def set_area_pixel(self, min_resol=1e-08, scheduler=None):
"""deprecated, obsolete"""
[docs]
@deprecated(details="This method is futile and will be removed in a future version."
" `Centroids.get_area_pixel` can be run without initialization.")
def set_area_approx(self, min_resol=1e-08):
"""deprecated, obsolete"""
[docs]
@deprecated(details="This method is futile and will be removed in a future version."
" `Centroids.get_dist_coast` can be run without initialization.")
def set_dist_coast(self, signed=False, precomputed=False, scheduler=None):
"""deprecated, obsolete"""
[docs]
@deprecated(details="This method has no effect and will be removed in a future version."
" In the current version of climada the geometry points of a `Centroids` object"
" cannot be removed as they are the backbone of the Centroids' GeoDataFrame.")
def empty_geometry_points(self):
""""deprecated, has no effect, which may be unexpected: no geometry points will be removed,
the centroids' GeoDataFrame is built on them!
"""
def _meta_to_lat_lon(meta):
"""Compute lat and lon of every pixel center from meta raster.
Parameters
----------
meta : dict
meta description of raster
Returns
-------
latitudes : np.ndarray
Latitudinal coordinates of pixel centers.
longitudes : np.ndarray
Longitudinal coordinates of pixel centers.
"""
xgrid, ygrid = u_coord.raster_to_meshgrid(meta['transform'], meta['width'], meta['height'])
return ygrid.ravel(), xgrid.ravel()
