Check if the centroid of a sub-basin is outside the polygon. If so create a new representative point #geoframe

check_geometry.py

#MIT License
#
#Copyright (c) 2017 Daniele Andreis
#Permission is hereby granted, free of charge, to any person obtaining a copy
#of this software and associated documentation files (the "Software"), to deal
#in the Software without restriction, including without limitation the rights
#to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#copies of the Software, and to permit persons to whom the Software is
#furnished to do so, subject to the following conditions:
#
#
#The above copyright notice and this permission notice shall be included in all
#copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
#SOFTWARE.

import rasterio as rasterio
from shapely.geometry import shape, mapping
import fiona
import os
from fiona.crs import from_epsg
from shapely import geometry, ops
import shutil


def cloneShapeFile(base_folder, source_shapefile_path, id):
    # Define the paths to the source (existing) shapefile and the destination (new) shapefile
    destination_shapefile_path = base_folder + "centroid_ID_old_" + id + ".shp"
    # Copy the source shapefile to the destination using shutil
    shutil.copy(source_shapefile_path, destination_shapefile_path)
    # Optional: If you want to rename the other associated files (.shx, .dbf, .prj, etc.), you can do so as well.
    shutil.copy(source_shapefile_path.replace('.shp', '.shx'), destination_shapefile_path.replace('.shp', '.shx'))
    shutil.copy(source_shapefile_path.replace('.shp', '.dbf'), destination_shapefile_path.replace('.shp', '.dbf'))
    shutil.copy(source_shapefile_path.replace('.shp', '.prj'), destination_shapefile_path.replace('.shp', '.prj'))

    print(f'Shapefile copied as {destination_shapefile_path}')


path = "net100/km1/geoframe/basins/"
doNewPoint = True
folders = []
for item in os.listdir(path):
    item_path = os.path.join(path, item)
    if os.path.isdir(item_path):
        folders.append(item)

print(folders)

for folder in folders:
    myCentroid = fiona.open(path + folder + "/centroid_ID_" + folder + ".shp")
    first = next(iter(myCentroid))
    point = shape(first['geometry'])
    polygons = fiona.open(path + folder + "/subbasins_complete_ID_" + folder + ".shp")
    first = next(iter(polygons))
    polygon = shape(first['geometry'])
    if not polygon.contains(point):
        cloneShapeFile(path + folder + "/", path + folder + "/centroid_ID_" + folder + ".shp", folder)
        print(folder)
        net_shp = path + folder + "/network_complete_ID_" + folder + ".shp"
        file_exists = os.path.exists(net_shp)
        # thanks to @Helios
        if file_exists and doNewPoint:
            net = fiona.open(net_shp)
            multiline = []
            for f in iter(net):
                multiline.append(shape(f['geometry']))
                multi_line = geometry.MultiLineString(multiline)
                merged_line = ops.linemerge(multi_line)
            new_centroid = merged_line.representative_point()
            output_shapefile = path + folder + "/centroid_ID_" + folder + ".shp"
            schema = {'geometry': 'Point',
              'properties': {'basinid': 'int', 'centrx': 'float', 'centry': 'float', 'elev_m': 'float',
                             'avgelev_m': 'float', 'area_km2': 'float', 'length_m': 'float',
                             'skyview': 'float'}}
            crs = from_epsg(32632)

            my_geometry = mapping(new_centroid)

            # Extract the x and y coordinates from the geometry
            x_coord, y_coord = my_geometry['coordinates']

            with rasterio.open(path + folder + "/dtm_"+folder+".asc") as src:
                # Use the geometry to get the pixel value(s) from the raster
                values = list(src.sample([my_geometry['coordinates']]))  # Returns a list of values at the geometry

                # If you want a single value for a point geometry, you can access it like this
                if my_geometry['type'] == 'Point':
                    elev = float(values[0][0])
                    print("Raster Values:", elev)

            with rasterio.open(path + folder + "/sky_"+folder+".asc") as src:
                # Use the geometry to get the pixel value(s) from the raster
                values = list(src.sample([my_geometry['coordinates']]))  # Returns a list of values at the geometry

                # If you want a single value for a point geometry, you can access it like this
                if my_geometry['type'] == 'Point':
                    sky = float(values[0][0])
                    print("Raster Values:", sky)


            with fiona.open(output_shapefile, 'w', 'ESRI Shapefile', schema, crs=crs) as output:
                feature = {
                    'geometry': my_geometry,
                    'properties': {'basinid': int(folder), 'centrx': float(x_coord), 'centry': float(y_coord) , 'elev_m': elev,
                                   'avgelev_m': first['properties']['avgelev_m'], 'area_km2': first['properties']['area_km2'],
                                   'length_m': first['properties']['length_m'],
                                   'skyview': sky}
                }
                output.write(feature)
        else:
            print("network not found")