Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
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()