LPT_functions.py

import dgpy
import geopandas as gpd
import numpy as np
from osgeo import ogr, osr, gdal
from shapely.geometry import Point


def initialize_zoi_map(path_mesh, mesh_proj, shp_file_name, field="fid"):
    bboxes = []
    groups = dgpy.dgGroupCollection(path_mesh)

    # Load reefmap and compute intersections between elements & reefs
    polys = dgpy.slimPolygonMeshIntersection(groups)

    ds = ogr.Open(shp_file_name, 0)
    lyr = ds.GetLayer()
    mesh_ref = ogr.osr.SpatialReference()
    mesh_ref.ImportFromProj4(mesh_proj)
    ctran = ogr.osr.CoordinateTransformation(lyr.GetSpatialRef(), mesh_ref)

    for i in range(lyr.GetFeatureCount()):
        # some weird convention difference between shp and gpkg happening ?
        if shp_file_name[-5:] == ".gpkg":
            i += 1

        f = lyr.GetFeature(i)
        ID_field = f.GetField(field)
        for ring in f.geometry():
            pts = ring.GetPoints()
            if not pts:
                continue
            ptsproj = np.array(
                np.asarray(ctran.TransformPoints(np.asarray(pts)))[:, :2]
            )
            bboxes.append(
                [
                    np.min(ptsproj[:, 0]),
                    np.min(ptsproj[:, 1]),
                    np.max(ptsproj[:, 0]),
                    np.max(ptsproj[:, 1]),
                ]
            )
            polys.addPolygon(ptsproj, ID_field)

    return polys


def random_points_in_polygon(number, polygon):
    """Generate random points in a given polygon"""
    points = []
    min_x, min_y, max_x, max_y = polygon.bounds
    i = 0
    while i < number:
        point = Point(np.random.uniform(min_x, max_x), np.random.uniform(min_y, max_y))
        if polygon.contains(point):
            points.append(point)
            i += 1
    return points  # returns list of shapely point


def gdf_seeding_points(n_pts, gdf_poly):
    """Return seeding points in a geodataframe"""
    gdf_points = gpd.GeoDataFrame()
    n_sites = len(gdf_poly)
    n_pts_per_poly = np.round(n_pts / n_sites)
    for ind in gdf_poly.index:
        points = random_points_in_polygon(n_pts_per_poly, gdf_poly.iloc[ind].geometry)
        id_col = np.repeat(gdf_poly.iloc[ind].ID, len(points))
        gdf_points_tmp = gpd.GeoDataFrame({"ID": id_col, "geometry": points})
        gdf_points = gdf_points.append(gdf_points_tmp, ignore_index=True)
    return gdf_points


def points_on_grid(x, y, nx, ny, minx, maxx, miny, maxy):
    """Locate a set of points on a regular grid"""
    mask_out = [all(tup) for tup in zip(x > minx, x < maxx, y > miny, y < maxy)]

    indx = np.int_(np.floor(((x[mask_out] - minx) / (maxx - minx)) * nx))
    indy = np.int_(np.floor(((y[mask_out] - miny) / (maxy - miny)) * ny))

    ind = np.ravel_multi_index(np.array([indx, indy]), (nx, ny), order="F")
    count = np.flip(np.reshape(np.bincount(ind, minlength=nx * ny), (ny, nx)), axis=0)

    return count


def export_tiff(rasterPath, minx, miny, pixelWidth, pixelHeight, array):
    """convert array to raster and export as a tiff file"""

    array = array[::-1]  # reverse array so the tif looks like the array

    rows = array.shape[0]
    cols = array.shape[1]

    driver = gdal.GetDriverByName("GTiff")
    outRaster = driver.Create(rasterPath, cols, rows, 1, gdal.GDT_Int32)
    outRaster.SetGeoTransform((minx, pixelWidth, 0, miny, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(32756)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()