Commit b8af7aa8 authored by Cecilia Nievas's avatar Cecilia Nievas
Browse files

Added function adjust_coord_of_polygons and its test

parent bb609e4c
......@@ -897,3 +897,71 @@ def retrieve_coords_to_adjust(cardinal_str):
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
......@@ -840,3 +840,58 @@ def test_determine_cardinal_point():
# 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
Polygon([(11.2, 39.2), (12.2, 39.2), (12.2, 40.2), (11.2, 40.2)]),
] # over
# Polygons in polygons2 are to the following cardinal points with respect to polygon1:
cardinal_pts = ["N", "S", "E", "W", "NE", "NW", "SE", "SW"]
# Results should be:
should_be_1 = [
Polygon([(lon_w, lat_s), (lon_e, lat_s), (lon_e, 41.1), (lon_w, 41.1)]),
Polygon([(lon_w, 40.1), (lon_e, 40.1), (lon_e, lat_n), (lon_w, lat_n)]),
Polygon([(lon_w, lat_s), (13.1, lat_s), (13.1, lat_n), (lon_w, lat_n)]),
Polygon([(12.1, lat_s), (lon_e, lat_s), (lon_e, lat_n), (12.1, lat_n)]),
Polygon([(lon_w, lat_s), (12.9, lat_s), (12.9, 40.9), (lon_w, 40.9)]),
Polygon([(11.9, lat_s), (lon_e, lat_s), (lon_e, 41.1), (11.9, 41.1)]),
Polygon([(lon_w, 39.9), (13.1, 39.9), (13.1, lat_n), (lon_w, lat_n)]),
Polygon([(12.1, 40.1), (lon_e, 40.1), (lon_e, lat_n), (12.1, lat_n)]),
]
should_be_2 = [
Polygon([(lon_w, 41.1), (lon_e, 41.1), (lon_e, 42.2), (lon_w, 42.2)]),
Polygon([(lon_w, 39.2), (lon_e, 39.2), (lon_e, 40.1), (lon_w, 40.1)]),
Polygon([(13.1, lat_s), (14.2, lat_s), (14.2, lat_n), (13.1, lat_n)]),
Polygon([(11.2, lat_s), (12.1, lat_s), (12.1, lat_n), (11.2, lat_n)]),
Polygon([(12.9, 40.9), (13.8, 40.9), (13.8, 41.8), (12.9, 41.8)]),
Polygon([(10.8, 41.1), (11.9, 41.1), (11.9, 42.2), (10.8, 42.2)]),
Polygon([(13.1, 38.8), (14.2, 38.8), (14.2, 39.9), (13.1, 39.9)]),
Polygon([(11.2, 39.2), (12.1, 39.2), (12.1, 40.1), (11.2, 40.1)]),
]
# Go one by one the elements of polygons2, calculate and assert:
for i in range(len(polygons2)):
function_poly1, function_poly2 = gdet_cr_ind.adjust_coord_of_polygons(
polygon1, polygons2[i], cardinal_pts[i]
)
assert function_poly1.bounds == should_be_1[i].bounds
assert function_poly1.area == should_be_1[i].area
assert function_poly2.bounds == should_be_2[i].bounds
assert function_poly2.area == should_be_2[i].area
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment