Source code for nes.methods.cell_measures

#!/usr/bin/env python

import numpy as np
from copy import deepcopy


[docs] def calculate_grid_area(self): """ Get coordinate bounds and call function to calculate the area of each cell of a grid. Parameters ---------- self : nes.Nes Source projection Nes Object. """ # Create bounds if they do not exist if self._lat_bnds is None or self._lon_bnds is None: self.create_spatial_bounds() # Get spatial number of vertices spatial_nv = self.lat_bnds['data'].shape[-1] # Reshape bounds if spatial_nv == 2: aux_shape = (self.lat_bnds['data'].shape[0], self.lon_bnds['data'].shape[0], 4) lon_bnds_aux = np.empty(aux_shape) lon_bnds_aux[:, :, 0] = self.lon_bnds['data'][np.newaxis, :, 0] lon_bnds_aux[:, :, 1] = self.lon_bnds['data'][np.newaxis, :, 1] lon_bnds_aux[:, :, 2] = self.lon_bnds['data'][np.newaxis, :, 1] lon_bnds_aux[:, :, 3] = self.lon_bnds['data'][np.newaxis, :, 0] lon_bnds = lon_bnds_aux del lon_bnds_aux lat_bnds_aux = np.empty(aux_shape) lat_bnds_aux[:, :, 0] = self.lat_bnds['data'][:, np.newaxis, 0] lat_bnds_aux[:, :, 1] = self.lat_bnds['data'][:, np.newaxis, 0] lat_bnds_aux[:, :, 2] = self.lat_bnds['data'][:, np.newaxis, 1] lat_bnds_aux[:, :, 3] = self.lat_bnds['data'][:, np.newaxis, 1] lat_bnds = lat_bnds_aux del lat_bnds_aux else: lon_bnds = self.lon_bnds['data'] lat_bnds = self.lat_bnds['data'] # Reshape bounds and assign as grid corner coordinates grid_corner_lon = deepcopy(lon_bnds).reshape(lon_bnds.shape[0]*lon_bnds.shape[1], lon_bnds.shape[2]) grid_corner_lat = deepcopy(lat_bnds).reshape(lat_bnds.shape[0]*lat_bnds.shape[1], lat_bnds.shape[2]) # Calculate cell areas grid_area = calculate_cell_area(grid_corner_lon, grid_corner_lat, earth_radius_minor_axis=self.earth_radius[0], earth_radius_major_axis=self.earth_radius[1]) return grid_area
[docs] def calculate_geometry_area(geometry_list, earth_radius_minor_axis=6356752.3142, earth_radius_major_axis=6378137.0): """ Get coordinate bounds and call function to calculate the area of each cell of a set of geometries. Parameters ---------- geometry_list : List List with polygon geometries. earth_radius_minor_axis : float Radius of the minor axis of the Earth. earth_radius_major_axis : float Radius of the major axis of the Earth. """ geometry_area = np.empty(shape=(len(geometry_list,))) for geom_ind in range(0, len(geometry_list)): # Calculate the area of each geometry in multipolygon and collection objects if geometry_list[geom_ind].type in ['MultiPolygon', 'GeometryCollection']: multi_geom_area = 0 for multi_geom_ind in range(0, len(geometry_list[geom_ind].geoms)): if geometry_list[geom_ind].geoms[multi_geom_ind].type == 'Point': continue geometry_corner_lon, geometry_corner_lat = geometry_list[geom_ind].geoms[multi_geom_ind].exterior.coords.xy geometry_corner_lon = np.array(geometry_corner_lon) geometry_corner_lat = np.array(geometry_corner_lat) geom_area = mod_huiliers_area(geometry_corner_lon, geometry_corner_lat) multi_geom_area += geom_area geometry_area[geom_ind] = multi_geom_area * earth_radius_minor_axis * earth_radius_major_axis # Calculate the area of each geometry else: geometry_corner_lon, geometry_corner_lat = geometry_list[geom_ind].exterior.coords.xy geometry_corner_lon = np.array(geometry_corner_lon) geometry_corner_lat = np.array(geometry_corner_lat) geom_area = mod_huiliers_area(geometry_corner_lon, geometry_corner_lat) geometry_area[geom_ind] = geom_area * earth_radius_minor_axis * earth_radius_major_axis return geometry_area
[docs] def calculate_cell_area(grid_corner_lon, grid_corner_lat, earth_radius_minor_axis=6356752.3142, earth_radius_major_axis=6378137.0): """ Calculate the area of each cell of a grid. Parameters ---------- grid_corner_lon : np.array Array with longitude bounds of grid. grid_corner_lat : np.array Array with longitude bounds of grid. earth_radius_minor_axis : float Radius of the minor axis of the Earth. earth_radius_major_axis : float Radius of the major axis of the Earth. """ # Calculate area for each grid cell n_cells = grid_corner_lon.shape[0] area = np.empty(shape=(n_cells,)) for i in range(0, n_cells): area[i] = mod_huiliers_area(grid_corner_lon[i], grid_corner_lat[i]) return area*earth_radius_minor_axis*earth_radius_major_axis
[docs] def mod_huiliers_area(cell_corner_lon, cell_corner_lat): """ Calculate the area of each cell according to Huilier's theorem. Reference: CDO (https://earth.bsc.es/gitlab/ces/cdo/) Parameters ---------- cell_corner_lon : np.array Longitude boundaries of each cell cell_corner_lat : np.array Latitude boundaries of each cell """ sum = 0 # Get points 0 (bottom left) and 1 (bottom right) in Earth coordinates point_0 = lon_lat_to_cartesian(cell_corner_lon[0], cell_corner_lat[0], earth_radius_major_axis=1) point_1 = lon_lat_to_cartesian(cell_corner_lon[1], cell_corner_lat[1], earth_radius_major_axis=1) point_0, point_1 = point_0[0], point_1[0] # Get number of vertices if cell_corner_lat[0] == cell_corner_lat[-1]: spatial_nv = len(cell_corner_lon) - 1 else: spatial_nv = len(cell_corner_lon) for i in range(2, spatial_nv): # Get point 2 (top right) in Earth coordinates point_2 = lon_lat_to_cartesian(cell_corner_lon[i], cell_corner_lat[i], earth_radius_major_axis=1) point_2 = point_2[0] # Calculate area of triangle between points 0, 1 and 2 sum += tri_area(point_0, point_1, point_2) # Copy to calculate area of next triangle if i == (spatial_nv - 1): point_1 = deepcopy(point_2) return sum
[docs] def tri_area(point_0, point_1, point_2): """ Calculate area between three points that form a triangle. Reference: CDO (https://earth.bsc.es/gitlab/ces/cdo/) Parameters ---------- point_0 : np.array Position of first point in cartesian coordinates. point_1 : np.array Position of second point in cartesian coordinates. point_2 : np.array Position of third point in cartesian coordinates. """ # Get length of side a (between point 0 and 1) tmp_vec = cross_product(point_0, point_1) sina = norm(tmp_vec) a = np.arcsin(sina) # Get length of side b (between point 0 and 2) tmp_vec = cross_product(point_0, point_2) sinb = norm(tmp_vec) b = np.arcsin(sinb) # Get length of side c (between point 1 and 2) tmp_vec = cross_product(point_2, point_1) sinc = norm(tmp_vec) c = np.arcsin(sinc) # Calculate area s = 0.5*(a+b+c) t = np.tan(s*0.5) * np.tan((s - a)*0.5) * np.tan((s - b)*0.5) * np.tan((s - c)*0.5) area = np.fabs(4.0 * np.arctan(np.sqrt(np.fabs(t)))) return area
[docs] def cross_product(a, b): """ Calculate cross product between two points. Parameters ---------- a : np.array Position of point a in cartesian coordinates. b : np.array Position of point b in cartesian coordinates. """ return [a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]]
[docs] def norm(cp): """ Normalize the result of the cross product operation. Parameters ---------- cp : np.array Cross product between two points. """ return np.sqrt(cp[0]*cp[0] + cp[1]*cp[1] + cp[2]*cp[2])
[docs] def lon_lat_to_cartesian(lon, lat, earth_radius_major_axis=6378137.0): """ Calculate lon, lat coordinates of a point on a sphere. Parameters ---------- lon : np.array Longitude values. lat : np.array Latitude values. earth_radius_major_axis : float Radius of the major axis of the Earth. """ lon_r = np.radians(lon) lat_r = np.radians(lat) x = earth_radius_major_axis * np.cos(lat_r) * np.cos(lon_r) y = earth_radius_major_axis * np.cos(lat_r) * np.sin(lon_r) z = earth_radius_major_axis * np.sin(lat_r) return np.column_stack([x, y, z])