Commit 996a5735 authored by Cecilia Nievas's avatar Cecilia Nievas
Browse files

Merge branch 'feature/ind04' into 'master'

Added subsequent functions to create cells around SERA industrial points

See merge request !14
parents 2574d4b2 3987a530
......@@ -654,3 +654,478 @@ def swell_cells_with_buffer(in_gdf, offset_lon, offset_lat):
out_gdf["geometry"] = new_polygons
return out_gdf
def define_threshold_angle(width_EW, width_NS, outangle="radians"):
"""This function defines the threshold angle to consider one cell at a particular cardinal
point with respect to another cell. The threshold angle depends only on the widths of the
cells, i.e. width_EW, width_NS. This threshold angle allows to distinguish, for example,
between a cell located to the north, east or north-east of another cell. Though the cell
widths are input in degrees, they are treated as Cartesian.
Args:
width_EW (float): Width of the cell in the east-west direction, in degrees, >0.
width_NS (float): Width of the cell in the north-south direction, in degrees, >0.
outangle (str): Indicates the units of the output angle, i.e. "radians" or "degrees".
Returns:
ang_thr (float): Threshold angle indicating:
- the angle between the north-south orientation and the limit lines to
consider a cell to be to the north or south of anoher cell (instead of
the north-east, or north-west, or south-east, or south-west)
- the angle between the east-west orientation and the limit lines to
consider a cell to be to the east or west of another cell (instead of
the north-east, or north-west, or south-east, or south-west)
It returns NaN if the input widths are not positive numbers.
"""
if np.logical_and(width_EW >= 0.0, width_NS >= 0.0):
ang_thr = np.arctan(width_NS / (2.0 * width_EW))
else:
print("ERROR!! width_EW AND width_NS NEED TO BE POSITIVE REAL NUMBERS")
ang_thr = np.nan
if outangle.lower() == "radians":
pass
elif outangle.lower() == "degrees":
ang_thr = np.rad2deg(ang_thr)
else:
print(
"ERROR!! outangle CAN ONLY BE radians OR degrees. "
"THRESHOLD ANGLE WILL BE RETURNED IN RADIANS"
)
return ang_thr
def angles_from_coords(p1x, p1y, p2x, p2y):
"""This function calculates the angle formed by an east-west line and a segment that starts
at point (p1x, p1y) and finishes at point (p2x, p2y).
Args:
p1x (array of floats) = X coordinates of Points labelled as 1.
p1y (array of floats) = Y coordinates of Points labelled as 1.
p2x (array of floats) = X coordinates of Points labelled as 2.
p2y (array of floats) = Y coordinates of Points labelled as 2.
Returns:
anglesfromcoord (array of floats): Angles formed by an E-W ("horizontal") line and a
segment that starts at P1 and finishes at P2, with
P1 being the intersection point between the two
lines. It will always be in the range [-pi, +pi].
If positive, it is measured counter-clockwise, if
negative, it is measured clockwise.
"""
anglesfromcoord = np.arctan2((p2y - p1y), (p2x - p1x))
return anglesfromcoord
def guarantee_radians(in_theta, theta_thr, inputangle):
"""This function converts in_theta and theta_thr into radians if inputangle is "degrees".
It returns the original in_theta and theta_thr if inputangle is "radians".
Args:
in_theta (float): Angles, in the units indicatd by inputangle
theta_thr (float): Angles, in the units indicatd by inputangle
Returns:
in_theta (float): Angles of in_theta, in radians
theta_thr (float): Angles of theta_thr, in radians
"""
if inputangle.lower() == "radians":
pass
elif inputangle.lower() == "degrees":
in_theta = np.deg2rad(in_theta)
theta_thr = np.deg2rad(theta_thr)
else:
print(
"ERROR!! outangle CAN ONLY BE radians OR degrees. "
"NOTHING DONE IN guarantee_radians!!!"
)
return in_theta, theta_thr
def determine_cardinal_point(in_theta, theta_thr, inputangle="radians"):
"""This function "translates" each angle in the array in_theta into a direction in terms of
"cardinal points". The array in_theta contains angles between an E-W ("horizontal") line
and another line, in the range [-pi, +pi], measured from the E-W ("horizontal") line. If
positive, it is measured counter-clockwise, if positive, it is measured clockwise. They are
expected to stem from the use of numpy.arctan2(), which yields an angle in this range,
under this same convention. The "translation" is done in terms of the tautological
definition of cardinal points (e.g. that 0 is east, pi/2 is north, pi and -pi are west,
-pi/2 is south) and an angle theta_thr that defines how much can the line deviate from,
for example, pi/2, and still be considered north. The exact ranges associated with each
of the eight possible outcomes (N, S, E, W, NE, NW, SE, SW) are defined within the code.
Args:
in_theta (array of floats): Angle between an E-W ("horizontal") line and another line,
in the range [-pi, +pi], measure from the E-W
("horizontal") line. If positive, it is measured counter-
clockwise, if negative, it is measured clockwise. It is
expected to stem from the use of numpy.arctan2(), which
yields an angle in this range, under this same convention.
If in_theta is outside the range [-pi, +pi], this function
will return an array with 'UU', meaning "unknown".
theta_thr (float): Angle used to divide all possible values of in_theta into east (E),
north-east (NE), north (N), north-west (NW), and so on. It needs to
be positive. The ranges that relate to each categories are defined
in the code.
inputangle (str): Units of in_theta and theta_thr, either "radians" or "degrees".
Returns:
out_card (array of str): Array of the same length as in_theta. Each individual value
can only be one of these values: N, S, E, W, NE, NW, SE, SW,
or U, the latter meaning "unknown", i.e. something is wrong
in the input.
"""
in_theta, theta_thr = guarantee_radians(in_theta, theta_thr, inputangle)
# Initialise output array with "unknown":
out_card = np.array(["UU" for i in range(len(in_theta))])
if theta_thr < 0.0: # theta_thr cannot be negative
return out_card
# East
out_card[np.where(np.logical_and(in_theta >= -theta_thr, in_theta <= theta_thr))[0]] = "E"
# North-east
out_card[
np.where(np.logical_and(in_theta > theta_thr, in_theta < (np.pi / 2.0 - theta_thr)))[0]
] = "NE"
# North
out_card[
np.where(
np.logical_and(
in_theta >= (np.pi / 2.0 - theta_thr), in_theta <= (np.pi / 2.0 + theta_thr)
)
)[0]
] = "N"
# North-west
out_card[
np.where(
np.logical_and(in_theta > (np.pi / 2.0 + theta_thr), in_theta < (np.pi - theta_thr))
)[0]
] = "NW"
# West (below, *1.00001 is the marging for floating point precision)
out_card[
np.where(
np.logical_or(
np.logical_and(in_theta >= (np.pi - theta_thr), in_theta <= np.pi * 1.00001),
np.logical_and(in_theta >= -np.pi * 1.00001, in_theta <= (-np.pi + theta_thr)),
)
)[0]
] = "W"
# South-west
out_card[
np.where(
np.logical_and(
in_theta > (-np.pi + theta_thr), in_theta < (-np.pi / 2.0 - theta_thr)
)
)[0]
] = "SW"
# South
out_card[
np.where(
np.logical_and(
in_theta >= (-np.pi / 2.0 - theta_thr), in_theta <= (-np.pi / 2.0 + theta_thr)
)
)[0]
] = "S"
# South-east
out_card[
np.where(np.logical_and(in_theta > (-np.pi / 2.0 + theta_thr), in_theta < -theta_thr))[
0
]
] = "SE"
return out_card
def retrieve_coords_to_adjust(cardinal_str):
"""This function returns the position of the coordinates that need to be changed in two
cells/tiles whose relative position is given by cardinal_str so that there is no overlap
or gap between them. The position refers to the west (0), south (1), east (2), north (3)
coordinates of shapely.polygon.bounds.
Args:
cardinal_str (str): N, S, E, W, NE, NW, SE, SW
Returns:
which_to_adj_1 (list): Position of the coordinate to be changed in the first cell
when the second cell is located to its cardinal_str. The
position refers to the west (0), south (1), east (2), north (3)
coordinates of shapely.polygon.bounds.
which_to_adj_2 (list): Position of the coordinate to be changed in the second cell
when it is located to the cardinal_str of the first cell.
The position refers to the west (0), south (1), east (2), north
(3) coordinates of shapely.polygon.bounds.
"""
if cardinal_str == "E":
which_to_adj_1 = [2] # East
which_to_adj_2 = [0] # West
elif cardinal_str == "NE":
which_to_adj_1 = [2, 3] # East & north
which_to_adj_2 = [0, 1] # West & south
elif cardinal_str == "N":
which_to_adj_1 = [3] # North
which_to_adj_2 = [1] # South
elif cardinal_str == "NW":
which_to_adj_1 = [3, 0] # North & west
which_to_adj_2 = [1, 2] # South & east
elif cardinal_str == "W":
which_to_adj_1 = [0] # West
which_to_adj_2 = [2] # East
elif cardinal_str == "SW":
which_to_adj_1 = [1, 0] # South & west
which_to_adj_2 = [3, 2] # North & east
elif cardinal_str == "S":
which_to_adj_1 = [1] # South
which_to_adj_2 = [3] # North
elif cardinal_str == "SE":
which_to_adj_1 = [1, 2] # South & east
which_to_adj_2 = [3, 0] # North & west
else:
which_to_adj_1 = [9] # error
which_to_adj_2 = [9] # error
return which_to_adj_1, which_to_adj_2
def adjust_coord_of_polygons(poly1, poly2, cardinal_str):
"""This function adjusts the geometry of the cells poly1 and poly2 so that there is no
overlap or gap between them. The relative position of poly2 with respect to poly1 is
given by cardinal_str, which can be N, S, E, W, NE, NW, SE, SW. This is used to determine
which coordinates to adjust, by calling retrieve_coords_to_adjust(). The new coordinate
is defined as the average of the corresponding coordinates of the two input polygons.
For example, if cardinal_str="NE", then poly2 is to the north-east of poly1, and the
coordinates that need to be adjusted are:
- the northern and eastern coordinates of poly1
- the southern and western coordinates of poly2
The new northern coordinate of poly1 and the new southern coordinate of poly2 will be the
average between the input northern coordinate of poly1 and the southern coordinate of poly2
Args:
poly1 (Shapely Polygon): First cell to adjust (cell implicitly a rectangle or square).
poly2 (Shapely Polygon): Second cell to adjust (cell implicitly a rectangle or square).
cardinal_str (str): Relative position of poly2 with respect to poly1. For example, if
cardinal_str="E", then poly2 is to the east of poly1, as determined
by the function determine_cardinal_point(). It can be N, S, E, W,
NE, NW, SE, SW.
Returns:
new_poly1 (Shapely Polygon): New geometry of poly1 so that there is no overlap or gap
between the two cells/tiles.
new_poly2 (Shapely Polygon): New geometry of poly2 so that there is no overlap or gap
between the two cells/tiles.
"""
# poly1_which and poly2_which will contain the positions of the coordinates to be changed
# (the positions refer to the order in which coords are displayed by polygon.bounds):
poly1_which, poly2_which = retrieve_coords_to_adjust(cardinal_str)
# Initialise the output polygons by copying the input ones:
new_poly1 = deepcopy(poly1)
new_poly2 = deepcopy(poly2)
# Go one by one the coordinates that need to be changed:
for j in range(len(poly1_which)):
# Old coordinate of poly1 to be changed:
coord1 = new_poly1.bounds[poly1_which[j]]
# Old coordinate of poly2 to be changed:
coord2 = new_poly2.bounds[poly2_which[j]]
# The new coordinate is an average of the two:
newcoord = (coord1 + coord2) / 2.0
newcoords1 = np.zeros([4])
newcoords2 = np.zeros([4])
for coordpos in range(4):
# Take all other coordinates from polygon.bounds:
if coordpos != poly1_which[j]:
newcoords1[coordpos] = new_poly1.bounds[coordpos]
else: # Take the new coordinate when the position needs to be changed
newcoords1[coordpos] = newcoord
# Take all other coordinates from polygon.bounds:
if coordpos != poly2_which[j]:
newcoords2[coordpos] = new_poly2.bounds[coordpos]
else: # Take the new coordinate when the position needs to be changed
newcoords2[coordpos] = newcoord
new_poly1 = define_cell_polygon(
newcoords1[0], newcoords1[1], newcoords1[2], newcoords1[3]
)
new_poly2 = define_cell_polygon(
newcoords2[0], newcoords2[1], newcoords2[2], newcoords2[3]
)
return new_poly1, new_poly2
def auto_adjust_overlaps_gaps(
in_gdf, overlaps_gdf, col_lon, col_lat, width_EW, width_NS, case, col_1="id_1", col_2="id_2"
):
"""This function resolves overlaps/gaps between cells by averaging the coordinates of the
corresponding boundaries. The overlaps existing in the cells of in_gdf are indicated in
overlaps_gdf, in which each row is an overlap case involving two cells. When case ==
"overlap", overlaps_gdf should contain the intersections of in_gdf itself, while when case
== "gap", overlaps_gdf should contain the intersections of an enlarged version of the cells
that allows to identify the potential gaps. When case == "overlap", overlaps_gdf is used
simply to identify neighbouring cells that overlap each other. When case == "gap",
overlaps_gdf is also used in combination with in_gdf to define whether there is a gap or
not, by subtracting the original cell geometries from in_gdf from the intersection of the
enlarged cells. If the resulting geometry is one Polygon, then there is a gap and this
function will fix it (by averaging the coordinates); if the resulting geometry is a
MultiPolygon then there is no gap to be fixed and the function moves on to analysing the
next case.
By calling angles_from_coords() and determine_cardinal_point(), the function identifies the
relative position of one cell with respect to the other (e.g. to the north, to the
south-east, etc.) to be able to decide which coordinates of the cells' geometries need to be
adjusted. The classification of the relative positions into cardinal points (eight cases,
north, north-east, east, south-east, south, south-west, west, north-west) is a function of
the widths of the cells (width_EW, width_NS), and is determined by calling
define_threshold_angle(). The adjustment of the cells' geometries is carried out by calling
adjust_coord_of_polygons().
Args:
in_gdf (Geopandas GeoDataFrame): GeoDataFrame with the cells geometries and original
centroids. Its "geometry" column contains the cells'
polygons. It also contains the columns "lon_w",
"lat_s", "lon_e", and "lat_n", with the four
coordinates that define each cell (floats). It is meant
to be the output of function
build_adjusted_cells_dataframe().
overlaps_gdf (Geopandas GeoDataFrame): GeoDataFrame with the geometry of the
intersections of the polygons in in_gdf (if case
== "overlap") or an expanded version of in_gdf
(if case == "gap"). It is meant to be the output
of function overlap_by_full_geom_intersection().
col_lon (str): Name of the column in the file of the aggregated exposure model that
contains longitudes.
col_lat (str): Name of the column in the file of the aggregated exposure model that
contains latitudes.
width_EW (float): Width of the cell in the East-West direction, in degrees, >0.
width_NS (float): Width of the cell in the North-south direction, in degrees, >0.
case (str): Indicates whether the function should deal with an "overlap" or "gap".
col_1 (str): Name of first column of overlaps_gdf to consider for comparing if the entry
already exists. Default: "id_1".
col_2 (str): Name of second column of overlaps_gdf to consider for comparing if the
entry already exists. Default: "id_2".
Returns:
out_gdf (Geopandas GeoDataFrame): GeoDataFrame with the same structure as was given as
input. The geometries of the cells for which
overlaps_gdf indicated there were overlaps are the
adjusted geometries calculated within this function.
The geometries of the rest of the cells are the same
as given in the input (no change).
gaps_found (bool): True if actual gaps were found and case == "gap". If case ==
"overlap", it will be False.
"""
# Check that case is either "overlap" or "gap":
if not (case == "overlap" or case == "gap"):
print("ERROR!! case can only be 'overlap' or 'gap' in auto_adjust_overlaps_gaps")
return {}
out_gdf = deepcopy(in_gdf)
# Calculate threshold angle that allows to define the relative position between overlapping
# cells:
ang_thr = define_threshold_angle(width_EW, width_NS)
# Calculate angles between te E-W line and a segment that starts at the centroid of the
# first cell and finishes at the centroid of the second cell, with the former centroid being
# the intersection between the two lines:
angles = angles_from_coords(
overlaps_gdf["%s_1" % (col_lon)].values,
overlaps_gdf["%s_1" % (col_lat)].values,
overlaps_gdf["%s_2" % (col_lon)].values,
overlaps_gdf["%s_2" % (col_lat)].values,
)
# Determine the cardinal points at which the second cell (col_2) is located with respect to
# the first cell (col_1):
cardinal_pts = determine_cardinal_point(angles, ang_thr, inputangle="radians")
# Go one by one the cases of overlaps/gaps found and adjust them:
gaps_found = False
for i in range(overlaps_gdf.shape[0]):
# IDs of the two cells involved:
id_1 = overlaps_gdf[col_1].values[i]
id_2 = overlaps_gdf[col_2].values[i]
if case == "gap":
how_many_intersections = max(
len(np.where(overlaps_gdf[col_1].values == id_1)[0]),
len(np.where(overlaps_gdf[col_1].values == id_2)[0]),
len(np.where(overlaps_gdf[col_2].values == id_1)[0]),
len(np.where(overlaps_gdf[col_2].values == id_2)[0]),
)
escape_condition = (how_many_intersections > 1) and (
(cardinal_pts[i] == "NW")
or (cardinal_pts[i] == "NE")
or (cardinal_pts[i] == "SW")
or (cardinal_pts[i] == "SE")
)
if escape_condition:
# If the cells are in diagonal with respect to one another and they participate
# of other cases of intersection, ignore:
continue
# Auxiliary GeoDataFrame with original cell geometries:
aux_orig = out_gdf[(out_gdf.id == id_1) | (out_gdf.id == id_2)]
# Auxiliary GeoDataFrame with just the geometry of this intersection:
aux_intersect = overlaps_gdf[
(overlaps_gdf.id_1 == id_1) & (overlaps_gdf.id_2 == id_2)
]
# Subtract the original cells geometry from the intersection:
subtraction = gpd.overlay(aux_intersect, aux_orig, how="difference")
if subtraction["geometry"].values[0].geom_type == "MultiPolygon":
# If MultiPolygon, no gap exists, thus no adjustment needed:
continue
gaps_found = True
# Identify the rows of these cells in out_gdf:
which_1 = np.where(out_gdf["id"].values == id_1)[0]
which_2 = np.where(out_gdf["id"].values == id_2)[0]
# Retrieve the original polygons of these cells:
polygon1 = out_gdf["geometry"].values[which_1[0]]
polygon2 = out_gdf["geometry"].values[which_2[0]]
# Calculate new polygons:
newpolygon1, newpolygon2 = adjust_coord_of_polygons(polygon1, polygon2, cardinal_pts[i])
# Write the new geometries to out_gdf:
out_gdf["geometry"].values[which_1[0]] = newpolygon1
out_gdf["geometry"].values[which_2[0]] = newpolygon2
# Write columns "lon_w", "lat_s", "lon_e", and "lat_n" of out_gdf, if they exist:
if (
"lon_w" in out_gdf.columns
and "lon_e" in out_gdf.columns
and "lat_s" in out_gdf.columns
and "lat_n" in out_gdf.columns
):
# Polygon 1:
out_gdf["lon_w"].values[which_1[0]] = newpolygon1.bounds[0]
out_gdf["lon_e"].values[which_1[0]] = newpolygon1.bounds[2]
out_gdf["lat_s"].values[which_1[0]] = newpolygon1.bounds[1]
out_gdf["lat_n"].values[which_1[0]] = newpolygon1.bounds[3]
# Polygon 2:
out_gdf["lon_w"].values[which_2[0]] = newpolygon2.bounds[0]
out_gdf["lon_e"].values[which_2[0]] = newpolygon2.bounds[2]
out_gdf["lat_s"].values[which_2[0]] = newpolygon2.bounds[1]
out_gdf["lat_n"].values[which_2[0]] = newpolygon2.bounds[3]
return out_gdf, gaps_found
......@@ -746,3 +746,486 @@ def test_swell_cells_with_buffer():
for i in range(len(expected_geometries)):
assert expected_geometries[i] == function_gdf["geometry"].values[i]
def test_determine_cardinal_point():
# Create input data for a realistic test:
theta_threshold = 25.0
angles_deg = np.array(
[
0.0,
24.0,
26.0,
64.0,
66.0,
90.0,
114.0,
116.0,
154.0,
156.0,
180.0,
185.0,
-0.1,
-24.0,
-26.0,
-64.0,
-66.0,
-90.0,
-114.0,
-116.0,
-154.0,
-156.0,
-180.0,
-185.0,
]
)
expected_output = np.array(
[
"E",
"E",
"NE",
"NE",
"N",
"N",
"N",
"NW",
"NW",
"W",
"W",
"UU",
"E",
"E",
"SE",
"SE",
"S",
"S",
"S",
"SW",
"SW",
"W",
"W",
"UU",
]
)
# Call function to test:
function_cardinal_pts = gdet_cr_ind.determine_cardinal_point(
angles_deg, theta_threshold, inputangle="degrees"
)
# Go one by one the elements of the output array:
for i in range(len(expected_output)):
assert function_cardinal_pts[i] == expected_output[i]
# Same test as above, but giving the angles in radians instead of degrees:
angles_rad = np.deg2rad(angles_deg)
theta_threshold_rad = np.deg2rad(theta_threshold)
# Call function to test:
function_cardinal_pts = gdet_cr_ind.determine_cardinal_point(
angles_rad, theta_threshold_rad, inputangle="radians"
)
# Go one by one the elements of the output array:
for i in range(len(expected_output)):
assert function_cardinal_pts[i] == expected_output[i]
# Test that all results are unknown ("UU") if the threshold angle is negative:
theta_threshold = -25.0
# Call function to test:
function_cardinal_pts = gdet_cr_ind.determine_cardinal_point(
angles_deg, theta_threshold, inputangle="degrees"
)
# Go one by one the elements of the output array:
for i in range(len(expected_output)):
assert function_cardinal_pts[i] == "UU"
def test_adjust_coord_of_polygons():
# Create datasets to test:
spacing = 1.0
# One polygon1, around which polygon2 will be "moved" to change their relative position:
lon_w = 12.0
lat_s = 40.0
lon_e = lon_w + spacing # 13.0
lat_n = lat_s + spacing # 41.0
polygon1 = Polygon([(lon_w, lat_s), (lon_e, lat_s), (lon_e, lat_n), (lon_w, lat_n)])
# An array of polygons at different positions with respect to polygon1:
polygons2 = [
Polygon([(lon_w, 41.2), (lon_e, 41.2), (lon_e, 42.2), (lon_w, 42.2)]), # gap
Polygon([(lon_w, 39.2), (lon_e, 39.2), (lon_e, 40.2), (lon_w, 40.2)]), # overlap
Polygon([(13.2, lat_s), (14.2, lat_s), (14.2, lat_n), (13.2, lat_n)]), # gap
Polygon([(11.2, lat_s), (12.2, lat_s), (12.2, lat_n), (11.2, lat_n)]), # overlap
Polygon([(12.8, 40.8), (13.8, 40.8), (13.8, 41.8), (12.8, 41.8)]), # overlap
Polygon([(10.8, 41.2), (11.8, 41.2), (11.8, 42.2), (10.8, 42.2)]), # gap
Polygon([(13.2, 38.8), (14.2, 38.8), (14.2, 39.8), (13.2, 39.8)]), # overlap