Skip to content
Snippets Groups Projects
Commit 59a0e44e authored by Sara Savanovic Djordjevic's avatar Sara Savanovic Djordjevic
Browse files

update: approach 4, grid of coordinates

parent 3b1ae240
No related branches found
No related tags found
1 merge request!6Clhp map
No preview for this file type
......@@ -13,7 +13,7 @@ def get_relation(self, body_of_water: str):
# Filter only polygons, exclude points and other feature types to reduce response size
polygon_data = geo_data[geo_data['geometry'].geom_type == 'Polygon']
if not polygon_data:
if polygon_data.empty:
raise ValueError("Failed to extract polygon data from file")
# Extract coordinates from polygons and create polygon objects
......@@ -26,22 +26,30 @@ def get_relation(self, body_of_water: str):
# Divide all polygons into sections
for polygon in polygons:
# Initialise grid
grid_lines = create_grid(polygon, cell_size=0.1)
hrz_lines = grid_lines[0] # Horizontal lines
vrt_lines = grid_lines[1] # Vertical lines
# Divide the length and with of polygon into a grid
grid_lines = create_grid_coords(polygon, cell_size=0.1)
hrz_lines = grid_lines[0] # Horizontal coordinates
vrt_lines = grid_lines[1] # Vertical coordinates
# Cut polygon into horizontal sections
for hrz_line in hrz_lines:
# Splits shape into upper and lower section by hrz_line
cut_poly_1 = cut_polygon_in_two(polygon, hrz_line)
# Split shape into upper and lower section as hrz_line as divider
cut_poly_1 = cut_polygon_in_two(polygon, hrz_line, True)
polygon_part = cut_poly_1[0] # Save upper horizontal section
if not polygon_part or not cut_poly_1[0]:
continue
# Cut each horizontal section into vertical sections
for vrt_line in vrt_lines:
cut_poly_2 = cut_polygon_in_two(polygon_part, vrt_line)
cut_poly_2 = cut_polygon_in_two(polygon_part, vrt_line, False)
if not cut_poly_2[0]:
continue
divided_map.extend(cut_poly_2[0]) # Append part to final list of shapes
# Set polygon_part to the remaining, un-split, horizontal section for next iteration
......@@ -70,66 +78,50 @@ def get_relation(self, body_of_water: str):
self.wfile.write(tiles_json.encode('utf-8'))
def cut_polygon_in_two(polygon, line):
# Check for invalid parameters
if not isinstance(polygon, Polygon) or not isinstance(line, LineString):
print("Inputs must be Shapely Polygon and LineString obects")
# Takes a polygon and divides its coordinates into two shapes, where divisor is a
# coordinate that defines the point of division
def cut_polygon_in_two(polygon: Polygon, divisor: float, horizontal: bool):
# Extract polygon exterior coordinates
exterior_coords = list(polygon.exterior.coords)
# Initialize lists to store points for the two shapes
shape_1 = []
shape_2 = []
# Initialize lists to store coordinates of new shapes after split
split_shape = []
remaining_shape = []
# Split points into separate shapes based on their location relative to the line
# Loop through points and check which side of the division line they are
for point in exterior_coords:
point = Point(point) # Convert coordinates to Shapely Point object
if line.distance(point) < 1e-10: # Floating point error tolerance
if line.contains(point): # Determine if point belongs to shape 1 or shape 2
shape_1.append(point)
if horizontal: # Horizontal split
if point.y < divisor:
split_shape.append(point)
else:
shape_2.append(point)
# Check if shapes are empty
if not shape_1 or not shape_2:
raise ValueError("One or both shapes are empty after cutting")
# Create Polygon objects from the lists of points
poly_1 = Polygon(shape_1)
poly_2 = Polygon(shape_2)
remaining_shape.append(point)
else: # Vertical split
if point.x < divisor:
split_shape.append(point)
else:
remaining_shape.append(point)
return poly_1, poly_2
# Check if polygons have enough coordinates
if len(split_shape) < 3 or len(remaining_shape) < 3:
print("Not enough coordinates to create valid polygons")
return None, None
# Append the first coordinate of the shapes to create closed loop
split_shape.append(split_shape[0]) # NB: may not be necessary?
remaining_shape.append(remaining_shape[0])
# From https://stackoverflow.com/questions/39338550/cut-a-polygon-with-two-lines-in-shapely
# Splits a polygon into two parts based on a dividing line
def cut_polygon_by_line(polygon, line):
if polygon.geom_type == 'Polygon':
polygon = [polygon]
merged = linemerge([poly.boundary for poly in polygon] + [line])
borders = unary_union(merged)
polygons = polygonize(borders)
return list(polygons)
return Polygon(split_shape), Polygon(remaining_shape) # Return split polygons as Shapely Polygon objects
# Create a grid of lines spaced by a cell size
def create_grid(polygon, cell_size):
# Generate grid of equally spaced x and y coordinates where the grid size is determined by cell_size
def create_grid_coords(polygon: Polygon, cell_size: float):
min_x, min_y, max_x, max_y = polygon.bounds
x_coords = np.arange(min_x, max_x, cell_size)
y_coords = np.arange(min_y, max_y, cell_size)
horizontal_lines = []
vertical_lines = []
# Vertical lines from top to bottom
for x in x_coords:
vertical_lines.append(LineString([(x, min_y), (x, max_y)]))
# Horizontal lines from top to bottom
for y in y_coords:
horizontal_lines.append(LineString([(min_x, y), (max_x, y)]))
if len(x_coords) == 0 or len(y_coords) == 0:
raise ValueError("No grid points generated")
if not horizontal_lines or not vertical_lines:
raise ValueError("List of horizontal or vertical lines is empty")
return x_coords, y_coords
return horizontal_lines, vertical_lines
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment