diff --git a/GDE_TOOLS_create_industrial_cells.py b/GDE_TOOLS_create_industrial_cells.py index b745a382c8c06e10aee71de3cc5dd486e135dbf7..5891b8495a4c7bad6e2c7bdc7fac5b4e9193bf99 100644 --- a/GDE_TOOLS_create_industrial_cells.py +++ b/GDE_TOOLS_create_industrial_cells.py @@ -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 diff --git a/tests/test_GDE_TOOLS_create_industrial_cells.py b/tests/test_GDE_TOOLS_create_industrial_cells.py index a99ea00c404b77bcda9c2a7909a9585a5da118ba..08b82c213961ea1b8cdd318d49c20b15205750fd 100644 --- a/tests/test_GDE_TOOLS_create_industrial_cells.py +++ b/tests/test_GDE_TOOLS_create_industrial_cells.py @@ -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 + 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 + + +def test_auto_adjust_overlaps_gaps(): + # Test that the function terminates if input variable "case" is not correct: + # ========================================================================== + assert {} == gdet_cr_ind.auto_adjust_overlaps_gaps([], [], "", "", 0.0, 0.0, "some_case") + + # General variables for all tests that follow: + col_lon = "LONGITUDE" + col_lat = "LATITUDE" + width_EW = 30.0 / (60.0 * 60.0) # 30 arcsec + width_NS = 30.0 / (60.0 * 60.0) # 30 arcsec + col_1 = "id_1" + col_2 = "id_2" + + # Test a case of overlaps: + # ======================== + # Define a GeoDataFrame with the cells' geometry: + smaller_spacing = 0.0082 + points_lon = np.array( + [ + 12.0, + 12.0, + 12.0 + smaller_spacing, + 12.0 + (smaller_spacing - 0.0001), + 12.0 + 2.0 * (smaller_spacing - 0.0001), + ] + ) + points_lat = np.array( + [ + 40.0, + 40.0 + smaller_spacing, + 40.0, + 40.0 + smaller_spacing, + 40.0 + 2.0 * (smaller_spacing - 0.0001), + ] + ) + longitudes_w = points_lon - width_EW / 2.0 + longitudes_e = points_lon + width_EW / 2.0 + latitudes_s = points_lat - width_NS / 2.0 + latitudes_n = points_lat + width_NS / 2.0 + geoms = [] + for i in range(0, len(points_lon)): + geoms.append( + Polygon( + [ + (longitudes_w[i], latitudes_s[i]), + (longitudes_e[i], latitudes_s[i]), + (longitudes_e[i], latitudes_n[i]), + (longitudes_w[i], latitudes_n[i]), + ] + ) + ) + + d = { + "id": np.array(["ID_%s" % (i) for i in range(1, len(geoms) + 1)]), + "LONGITUDE": points_lon, + "LATITUDE": points_lat, + "lon_w": longitudes_w, + "lon_e": longitudes_e, + "lat_s": latitudes_s, + "lat_n": latitudes_n, + "geometry": geoms, + } + cells_gdf = gpd.GeoDataFrame(d, geometry="geometry", crs="EPSG:4326") + + # Define a GeoDataFrame with the cells' intersections: + ids = {} + ids[1] = ["ID_1", "ID_1", "ID_1", "ID_2", "ID_2", "ID_3", "ID_4"] + ids[2] = ["ID_2", "ID_3", "ID_4", "ID_3", "ID_4", "ID_4", "ID_5"] + geoms = [ + Polygon( + [ + (12.0 - width_EW / 2.0, 40.0 + smaller_spacing - width_NS / 2.0), # IDs 1 & 2 + (12.0 + width_EW / 2.0, 40.0 + smaller_spacing - width_NS / 2.0), + (12.0 + width_EW / 2.0, 40.0 + width_NS / 2.0), + (12.0 - width_EW / 2.0, 40.0 + width_NS / 2.0), + ] + ), + Polygon( + [ + (12.0 + smaller_spacing - width_EW / 2.0, 40.0 - width_NS / 2.0), # IDs 1 & 3 + (12.0 + width_EW / 2.0, 40.0 - width_NS / 2.0), + (12.0 + width_EW / 2.0, 40.0 + width_NS / 2.0), + (12.0 + smaller_spacing - width_EW / 2.0, 40.0 + width_NS / 2.0), + ] + ), + Polygon( + [ + ( + 12.0 + (smaller_spacing - 0.0001) - width_EW / 2.0, + 40.0 + smaller_spacing - width_NS / 2.0, + ), # IDs 1 & 4 + (12.0 + width_EW / 2.0, 40.0 + smaller_spacing - width_NS / 2.0), + (12.0 + width_EW / 2.0, 40.0 + width_NS / 2.0), + (12.0 + (smaller_spacing - 0.0001) - width_EW / 2.0, 40.0 + width_NS / 2.0), + ] + ), + Polygon( + [ + ( + 12.0 + smaller_spacing - width_EW / 2.0, + 40.0 + smaller_spacing - width_NS / 2.0, + ), # IDs 2 & 3 + (12.0 + width_EW / 2.0, 40.0 + smaller_spacing - width_NS / 2.0), + (12.0 + width_EW / 2.0, 40.0 + width_NS / 2.0), + (12.0 + smaller_spacing - width_EW / 2.0, 40.0 + width_NS / 2.0), + ] + ), + Polygon( + [ + ( + 12.0 + (smaller_spacing - 0.0001) - width_EW / 2.0, + 40.0 + smaller_spacing - width_NS / 2.0, + ), # IDs 2 & 4 + (12.0 + width_EW / 2.0, 40.0 + smaller_spacing - width_NS / 2.0), + (12.0 + width_EW / 2.0, 40.0 + smaller_spacing + width_NS / 2.0), + ( + 12.0 + (smaller_spacing - 0.0001) - width_EW / 2.0, + 40.0 + smaller_spacing + width_NS / 2.0, + ), + ] + ), + Polygon( + [ + ( + 12.0 + smaller_spacing - width_EW / 2.0, + 40.0 + smaller_spacing - width_NS / 2.0, + ), # IDs 3 & 4 + ( + 12.0 + (smaller_spacing - 0.0001) + width_EW / 2.0, + 40.0 + smaller_spacing - width_NS / 2.0, + ), + (12.0 + (smaller_spacing - 0.0001) + width_EW / 2.0, 40.0 + width_NS / 2.0), + (12.0 + smaller_spacing - width_EW / 2.0, 40.0 + width_NS / 2.0), + ] + ), + Polygon( + [ + ( + 12.0 + 2.0 * (smaller_spacing - 0.0001) - width_EW / 2.0, + 40.0 + 2.0 * (smaller_spacing - 0.0001) - width_NS / 2.0, + ), # IDs 4 & 5 + ( + 12.0 + (smaller_spacing - 0.0001) + width_EW / 2.0, + 40.0 + 2.0 * (smaller_spacing - 0.0001) - width_NS / 2.0, + ), + ( + 12.0 + (smaller_spacing - 0.0001) + width_EW / 2.0, + 40.0 + smaller_spacing + width_NS / 2.0, + ), + ( + 12.0 + 2.0 * (smaller_spacing - 0.0001) - width_EW / 2.0, + 40.0 + smaller_spacing + width_NS / 2.0, + ), + ] + ), + ] + + cols = ["LONGITUDE", "LATITUDE", "lon_w", "lon_e", "lat_s", "lat_n"] + aux_dict = {} + + for col_name in cols: + for val in range(1, 3): + aux_list = [] + for i in range(0, len(ids[val])): + aux_list.append(cells_gdf[cells_gdf.id == ids[val][i]][col_name].values[0]) + aux_dict["%s_%s" % (col_name, str(val))] = np.array(aux_list) + d = {"id_1": ids[1], "id_2": ids[2]} + for key in aux_dict.keys(): + d[key] = aux_dict[key] + d["geometry"] = geoms + intersections_gdf = gpd.GeoDataFrame(d, geometry="geometry", crs="EPSG:4326") + + # Expected output: + expected_out_gdf = deepcopy(cells_gdf) + # Manual adustment done by row of the intersections GeoDataFrame (i.e. the order matters): + # IDs 1 & 2 + expected_out_gdf["lat_n"].values[0] = ( + expected_out_gdf["lat_n"].values[0] + expected_out_gdf["lat_s"].values[1] + ) / 2.0 + expected_out_gdf["lat_s"].values[1] = expected_out_gdf["lat_n"].values[0] + # IDs 1 & 3: + expected_out_gdf["lon_e"].values[0] = ( + expected_out_gdf["lon_e"].values[0] + expected_out_gdf["lon_w"].values[2] + ) / 2.0 + expected_out_gdf["lon_w"].values[2] = expected_out_gdf["lon_e"].values[0] + # IDs 1 & 4: + expected_out_gdf["lat_n"].values[0] = ( + expected_out_gdf["lat_n"].values[0] + expected_out_gdf["lat_s"].values[3] + ) / 2.0 + expected_out_gdf["lat_s"].values[3] = expected_out_gdf["lat_n"].values[0] + expected_out_gdf["lon_e"].values[0] = ( + expected_out_gdf["lon_e"].values[0] + expected_out_gdf["lon_w"].values[3] + ) / 2.0 + expected_out_gdf["lon_w"].values[3] = expected_out_gdf["lon_e"].values[0] + # IDs 2 & 3: + expected_out_gdf["lon_e"].values[1] = ( + expected_out_gdf["lon_e"].values[1] + expected_out_gdf["lon_w"].values[2] + ) / 2.0 + expected_out_gdf["lon_w"].values[2] = expected_out_gdf["lon_e"].values[1] + expected_out_gdf["lat_s"].values[1] = ( + expected_out_gdf["lat_s"].values[1] + expected_out_gdf["lat_n"].values[2] + ) / 2.0 + expected_out_gdf["lat_n"].values[2] = expected_out_gdf["lat_s"].values[1] + # IDs 2 & 4: + expected_out_gdf["lon_e"].values[1] = ( + expected_out_gdf["lon_e"].values[1] + expected_out_gdf["lon_w"].values[3] + ) / 2.0 + expected_out_gdf["lon_w"].values[3] = expected_out_gdf["lon_e"].values[1] + # IDs 3 & 4: + expected_out_gdf["lat_n"].values[2] = ( + expected_out_gdf["lat_n"].values[2] + expected_out_gdf["lat_s"].values[3] + ) / 2.0 + expected_out_gdf["lat_s"].values[3] = expected_out_gdf["lat_n"].values[2] + # IDs 4 & 5: + expected_out_gdf["lat_n"].values[3] = ( + expected_out_gdf["lat_n"].values[3] + expected_out_gdf["lat_s"].values[4] + ) / 2.0 + expected_out_gdf["lat_s"].values[4] = expected_out_gdf["lat_n"].values[3] + expected_out_gdf["lon_e"].values[3] = ( + expected_out_gdf["lon_e"].values[3] + expected_out_gdf["lon_w"].values[4] + ) / 2.0 + expected_out_gdf["lon_w"].values[4] = expected_out_gdf["lon_e"].values[3] + + new_geoms = [] + for i in range(0, expected_out_gdf.shape[0]): + lon_w = expected_out_gdf["lon_w"].values[i] + lon_e = expected_out_gdf["lon_e"].values[i] + lat_s = expected_out_gdf["lat_s"].values[i] + lat_n = expected_out_gdf["lat_n"].values[i] + new_geoms.append( + Polygon([(lon_w, lat_s), (lon_e, lat_s), (lon_e, lat_n), (lon_w, lat_n)]) + ) + expected_out_gdf["geometry"] = new_geoms + + # Output from function to be tested: + function_out_gdf, flag = gdet_cr_ind.auto_adjust_overlaps_gaps( + cells_gdf, + intersections_gdf, + col_lon, + col_lat, + width_EW, + width_NS, + "overlap", + col_1=col_1, + col_2=col_2, + ) + assert flag == False + + for rownum in range(0, expected_out_gdf.shape[0]): + bounds_expected = expected_out_gdf["geometry"].values[rownum].bounds + bounds_function = function_out_gdf["geometry"].values[rownum].bounds + for position in range(0, 4): + assert round(bounds_expected[position], 12) == round(bounds_function[position], 12) + np.testing.assert_allclose( + expected_out_gdf["lon_w"].values, function_out_gdf["lon_w"].values, rtol=0.0, atol=1e-08 + ) + np.testing.assert_allclose( + expected_out_gdf["lon_e"].values, function_out_gdf["lon_e"].values, rtol=0.0, atol=1e-08 + ) + np.testing.assert_allclose( + expected_out_gdf["lat_s"].values, function_out_gdf["lat_s"].values, rtol=0.0, atol=1e-08 + ) + np.testing.assert_allclose( + expected_out_gdf["lat_n"].values, function_out_gdf["lat_n"].values, rtol=0.0, atol=1e-08 + ) + + # Test a case of gaps: + # ==================== + # The output GeoDataFrame from the test of overlaps (expected_out_gdf) contains gaps and + # is used herein to test case="gap" (this would be the normal flow of the main code too): + # The new input is: + cells_gdf = deepcopy(expected_out_gdf) + cells_offset_gdf = gdet_cr_ind.swell_cells_with_buffer( + cells_gdf, 0.25 * width_EW, 0.25 * width_NS + ) + intersections_gdf = gdet_cr_ind.overlap_by_full_geom_intersection( + cells_offset_gdf, "id_1", "id_2" + ) + # The new expected output is: + # IDs 1 & 2 + expected_out_gdf["lat_n"].values[0] = ( + expected_out_gdf["lat_n"].values[0] + expected_out_gdf["lat_s"].values[1] + ) / 2.0 + expected_out_gdf["lat_s"].values[1] = expected_out_gdf["lat_n"].values[0] + # IDs 1 & 3: + expected_out_gdf["lon_e"].values[0] = ( + expected_out_gdf["lon_e"].values[0] + expected_out_gdf["lon_w"].values[2] + ) / 2.0 + expected_out_gdf["lon_w"].values[2] = expected_out_gdf["lon_e"].values[0] + + new_geoms = [] + for i in range(0, expected_out_gdf.shape[0]): + lon_w = expected_out_gdf["lon_w"].values[i] + lon_e = expected_out_gdf["lon_e"].values[i] + lat_s = expected_out_gdf["lat_s"].values[i] + lat_n = expected_out_gdf["lat_n"].values[i] + new_geoms.append( + Polygon([(lon_w, lat_s), (lon_e, lat_s), (lon_e, lat_n), (lon_w, lat_n)]) + ) + expected_out_gdf["geometry"] = new_geoms + + # Output from function to be tested: + function_out_gdf, flag = gdet_cr_ind.auto_adjust_overlaps_gaps( + cells_gdf, + intersections_gdf, + col_lon, + col_lat, + width_EW, + width_NS, + "gap", + col_1=col_1, + col_2=col_2, + ) + assert flag == True + + for rownum in range(0, expected_out_gdf.shape[0]): + bounds_expected = expected_out_gdf["geometry"].values[rownum].bounds + bounds_function = function_out_gdf["geometry"].values[rownum].bounds + for position in range(0, 4): + assert round(bounds_expected[position], 12) == round(bounds_function[position], 12) + np.testing.assert_allclose( + expected_out_gdf["lon_w"].values, function_out_gdf["lon_w"].values, rtol=0.0, atol=1e-08 + ) + np.testing.assert_allclose( + expected_out_gdf["lon_e"].values, function_out_gdf["lon_e"].values, rtol=0.0, atol=1e-08 + ) + np.testing.assert_allclose( + expected_out_gdf["lat_s"].values, function_out_gdf["lat_s"].values, rtol=0.0, atol=1e-08 + ) + np.testing.assert_allclose( + expected_out_gdf["lat_n"].values, function_out_gdf["lat_n"].values, rtol=0.0, atol=1e-08 + )