Cartography/prettymaps
1
0
Fork 0
mirror of https://github.com/marceloprates/prettymaps.git synced 2025-08-01 12:20:37 +02:00
Code Issues Projects Releases Wiki Activity

Implement Plotting coastline with just perimeter and change CRS for loaded polygons

Browse Source 
Now using WGS84 shape file https://osmdata.openstreetmap.de/download/water-polygons-split-4326.zip
This commit is contained in:
Colm McDonald
2021-08-31 20:54:38 +01:00
parent ab997b48d5
commit 577bd33393
1 changed files with 8 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

17
prettymaps/fetch.py
View File

@@ -1,6 +1,7 @@
from functools import reduce
import osmnx as ox
from osmnx import utils_geo
import numpy as np
from shapely.geometry import Point, Polygon, MultiPolygon, MultiLineString
from shapely.ops import unary_union
@@ -30,21 +31,19 @@ def get_perimeter(query, by_osmid = False, **kwargs):
def get_coast(perimeter = None, point = None, radius = None, tags = {}, perimeter_tolerance = 0, union = True, circle = True, dilate = 0, file_location = None):
if perimeter is not None:
# Boundary defined by polygon (perimeter)
geometries = ox.geometries_from_polygon(
unary_union(perimeter.geometry).buffer(perimeter_tolerance) if perimeter_tolerance > 0 else unary_union(perimeter.geometry),
tags = {tags: True} if type(tags) == str else tags
)
bbox=perimeter.to_crs(3174)
bbox=bbox.buffer(perimeter_tolerance+dilate)
bbox=bbox.to_crs(4326)
bbox=bbox.envelope
geometries = read_file(file_location, bbox=bbox)
perimeter = unary_union(ox.project_gdf(perimeter).geometry)
elif (point is not None) and (radius is not None):
# Boundary defined by circle with radius 'radius' around point
bbox=GeoDataFrame(geometry = [Point(point[::-1])], crs = 4326)
bbox=bbox.to_crs(3174)
bbox=bbox.buffer(radius+dilate)
bbox=bbox.envelope
north,south,west,east=utils_geo.bbox_from_point(point, dist=radius+dilate)
bbox=(west,south,east,north)
# Load the polygons for the coastline from a file
geometries=read_file(file_location, bbox=bbox)
geometries=geometries.to_crs("epsg:4326")
perimeter = get_boundary(point, radius, geometries.crs, circle = circle, dilate = dilate)
# Project GDF