Create additional layers in postprocessing step (The sea is not queryable in OSM)

[samdabadei]

It would be nice if you could create additional layers from existing gdfs-entries with the postprocessing function.
I imagine it in a way that you can give the layers a style but the layers only gets added in the postprocessing step. I tried to give an empty layer with in the layers-dict but this yields an error and without the layer beeing added in the layers-dict it wont be rendered.

Why do I need it? Apparently there is no tag for the sea, ocean in OSM (appart bays, etc). I think seas are therefore created implicitly with the cutout of the costlines. But i want to give the whole ocean a special styling. (For example: Try to render the island Lehtisaari in Helsinki. You will see that the ocean is rendered as if it would be the background.)

Thanks for tips to solve this without the mentioned feature :)

[samdabadei]

To just fill the costline unfortunately is also not possible. Im working with this:

darkstyle_helsinke = {layers: {perimeters:{}, costline:{tags:costline:True}}}, 
style:{ "perimeter": {"fill": True, "zorder": -1,  "fc": "#a8e1e6", "ec": "#2F3737",  "hatch_c": "#9bc3d4",  "hatch": "ooo...",  "lw": 2,},
        "costline": {"fill":True, "fc": "#9E5B20", "ec": "#1B3245", "lw": 1, "zorder": 0},}}
plot = pm.plot(
    "Lehtisaari, Helsinki",
    **darkstyle_helsinki,
    radius=500,
    credit=False,
    preset=None
)

[samdabadei]

if someone is interested in a dirty postprocessing-functions to close all costlines to polygones so that you can differentiate between land and ocean here you go:

`
import geopandas
from shapely import geometry, ops
from math import radians, cos, sin, asin, sqrt

def costline_processing(gdf):

def distance(lat1, lat2, lon1, lon2):
# The math module contains a function named
# radians which converts from degrees to radians.
lon1 = radians(lon1)
lon2 = radians(lon2)
lat1 = radians(lat1)
lat2 = radians(lat2)

# Haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2

c = 2 * asin(sqrt(a))

# Radius of earth in kilometers. Use 3956 for miles
r = 6371

# calculate the result
return (c * r)

def distance_two_points(p1, p2):
return distance(p1.xy[0], p2.xy[0], p1.xy[1], p2.xy[1]) * 1000

def create_linestring_from_xy_tuples(list_of_tuples):
    point_list = []

    for el in list_of_tuples:
        for x, y in zip(*el):
            point_list.append([x,y])

    return geometry.LineString(point_list)

def get_point_of_line(line, index):
    return (line.xy[0][index], line.xy[1][index])


frame = geometry.polygon.orient(gdf["perimeter"].geometry[0], sign=1)
frame_line_segments = geometry.MultiLineString([geometry.LineString([p1, p2]) for p1, p2 in zip(zip(frame.exterior.xy[0][:-1],frame.exterior.xy[1][:-1]),
                                                                                                zip(frame.exterior.xy[0][1:],frame.exterior.xy[1][1:]))])
geoms = gdf["costline"].geometry.copy()

#add geometries of costline to resprective list
geoms_poly = []
geoms_lines = []

for el in geoms:
    if type(el) == geometry.polygon.Polygon:
        geoms_poly.append(el)
    elif type(el) == geometry.linestring.LineString:
        geoms_lines.append(el)
    elif type(el) == geometry.multilinestring.MultiLineString:
        geoms_lines.extend([ell for ell in el])
    else:
        raise TypeError(f"Unknown geometry type: {type(el)}")



#merge lines if they are connected
merged_lines = ops.linemerge(geometry.MultiLineString(geoms_lines))

#merged_lines = geometry.MultiLineString([el.intersection(frame) for el in merged_lines])

#check for poligons that are closed
indexes_to_remove = []
for ind, el in enumerate(merged_lines):
    if get_point_of_line(el, 0) == get_point_of_line(el, -1):
        geoms_poly.append(geometry.Polygon(el))
        indexes_to_remove.append((ind))

merged_lines = geometry.MultiLineString([el for ind, el in enumerate(merged_lines) if ind not in indexes_to_remove])

# close lines for poligons over frame
closed_lines = []

for idx, el in enumerate(merged_lines):
    start = geometry.Point(el.xy[0][0], el.xy[1][0])
    end = geometry.Point(el.xy[0][-1], el.xy[1][-1])

    distances_start = [ops.nearest_points(line_i, start)[0].distance(start) for line_i in frame_line_segments]
    line_el_start = distances_start.index(min(distances_start))
    line_point_start = ops.nearest_points(frame_line_segments[line_el_start], start)[0]

    distances_end = [ops.nearest_points(line_i, end)[0].distance(end) for line_i in frame_line_segments]
    line_el_end = distances_end.index(min(distances_end))
    line_point_end = ops.nearest_points(frame_line_segments[line_el_end], end)[0]

    index_order = list(range(line_el_end, 4)) + list(range(0, line_el_end))

    point = line_point_end
    new_line = el
    for line_index in index_order:

        line_comp = geometry.LineString([[point.xy[0][0], point.xy[1][0]],
                                         [frame_line_segments[line_index].xy[0][-1], frame_line_segments[line_index].xy[1][-1]]])
        if line_comp.distance(line_point_start)<10e-10:
            new_line = create_linestring_from_xy_tuples([new_line.xy, start.xy])
            break
        else:
            point = geometry.Point(frame_line_segments[line_index].xy[0][-1], frame_line_segments[line_index].xy[1][-1])
            new_line = create_linestring_from_xy_tuples([new_line.xy, point.xy])

    if get_point_of_line(new_line, 0) != get_point_of_line(new_line, -1):
        new_line = create_linestring_from_xy_tuples([new_line.xy, start.xy])

    closed_lines.append(new_line)

for line in closed_lines:
    geoms_poly.append(geometry.Polygon(line))

#check if some poligons are sea-poligons
sea_idx = []
for idx, poly in enumerate(geoms_poly):
    share_area = [poly.intersection(poly_i).area for poly_i in geoms_poly]
    if sum([1 if el > 0 else 0 for el in share_area]) == len(geoms_poly):
        sea_idx.append(idx)

if len(sea_idx) > 0:
    #all polies are sea-poligons ore within the only land poligon
    if len(sea_idx) == len(geoms_poly):
        land_poly = frame

        for poly_idx in sea_idx:
            land_poly = land_poly.intersection(geoms_poly[poly_idx])

        geoms_poly = [land_poly]

    else:
        land_idx = [idx for idx in range(len(geoms_poly)) if idx not in sea_idx]

        #check in wich poligon the sea-poligon is

        respective_land_poly_idx = []
        factor = 0.9999
        for sea_poly_idx in sea_idx:
            for land_poly_idx in land_idx:
                if geoms_poly[sea_poly_idx].intersection(geoms_poly[land_poly_idx]).area < factor*geoms_poly[land_poly_idx].area:
                    respective_land_poly_idx.append(land_poly_idx)

        #add land poligons with no sea-poligons inside
        geoms_poly_new = [geoms_poly[idx] for idx in land_idx if idx not in respective_land_poly_idx]

        cut_land_dict = {}
        for idx, res_land_poly_idx in enumerate(respective_land_poly_idx):
            land = cut_land_dict.get(f"{res_land_poly_idx}", geoms_poly[res_land_poly_idx])
            land = land.intersection(geoms_poly[sea_idx[idx]])
            cut_land_dict[f"{res_land_poly_idx}"] = land

        for key, el in cut_land_dict.items():
            geoms_poly_new.append(el)

        geoms_poly = geoms_poly_new

d = {"geometry": geoms_poly}
gdf["costline"] = geopandas.GeoDataFrame(d, crs=gdf["perimeter"].crs)


return gdf

`