GISfunctions.py 17 KB
Newer Older
luroth's avatar
luroth committed
1
2
3
4
5
6
7
8
9
10
##############################################################################
#
# Common functions
#
# Author: Lukas Roth (lukas.roth@usys.ethz.ch)
#
##############################################################################

# Dependency
import math
luroth's avatar
luroth committed
11

luroth's avatar
luroth committed
12
13
14
15
import numpy as np
import stl
from osgeo import gdal
from shapely.geometry import Point, Polygon
luroth's avatar
luroth committed
16
17
import gdal
import pandas as pd
luroth's avatar
luroth committed
18
import geopandas as gpd
luroth's avatar
luroth committed
19
20
import geojson
from matplotlib import path
luroth's avatar
luroth committed
21
22
from scipy.spatial import Voronoi
from sklearn.cluster import AgglomerativeClustering
luroth's avatar
luroth committed
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

def ray_plane_intersect(ray_point, ray_vector, plane_point, plane_normal, epsilon=1e-64):
    """Intersect ray with plane

    :param ray_point: Numpy array with coordinates of point on ray
    :param ray_vector: Numpy array with direction of ray
    :param plane_point: Numpy array with coordinates of point on plane
    :param plane_normal: Numpy array with normal of plane
    :param epsilon: Epsilon value as float. Intersection distances with a lower value than epsilon are ignored
    :return:
    """

    ndotu = plane_normal.dot(ray_vector)

    if abs(ndotu) < epsilon:
        # no intersection or line is within plane
        return None
    else:
        # Trace ray. Based on:
        # Hughes, J. F. (2014). Computer graphics : principles and practice. Upper Saddle River, N.J: Addison-Wesley, 3rd ed. edition.
        w = ray_point - plane_point
        si = -plane_normal.dot(w) / ndotu
        psi = w + si * ray_vector + plane_point
        return psi


luroth's avatar
luroth committed
49
def image_rgb_cut(point, polygon):
luroth's avatar
luroth committed
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
104
105
106
107
108
109
110
111
112
113
114
115
116
    """Test if point is inside poylgon. Using

    :param point: Tuple with coordinates of point to test
    :param polygon: List of tuples with coordinates of polygon to test
    :return:
    """
    point = Point(point[0], point[1])
    polygon = Polygon(polygon)
    return polygon.contains(point)


def find_bounding_box_mesh(mesh):
    """Find STL mesh bounding box in xyz

    :param mesh: STL mesh
    :return: minx, maxx, miny, maxy, minz, maxz
    """

    minx = maxx = miny = maxy = minz = maxz = None
    for p in mesh.points:
        maxx = max(value for value in [p[stl.Dimension.X], maxx] if value is not None)
        minx = min(value for value in [p[stl.Dimension.X], minx] if value is not None)
        maxy = max(value for value in [p[stl.Dimension.Y], maxy] if value is not None)
        miny = min(value for value in [p[stl.Dimension.Y], miny] if value is not None)
        maxz = max(value for value in [p[stl.Dimension.Z], maxz] if value is not None)
        minz = min(value for value in [p[stl.Dimension.Z], minz] if value is not None)

    return minx, maxx, miny, maxy, minz, maxz


def get_triangle_normale(points):
    """Calculates normale of a plane defined by three edge points

    :param points: List of 3 points
    :return: Normale
    """

    a = np.array(points[0])
    b = np.array(points[1])
    c = np.array(points[2])

    v1 = a - b
    v2 = a - c
    normale = np.cross(v1, v2)

    return normale


def rotation_matrix(axis, theta):
    """ Create 3-D rotation matrix
    Based on Euler–Rodrigues formula (https://en.wikipedia.org/wiki/Euler%E2%80%93Rodrigues_formula)

    :param axis: Axis to rotate on
    :param theta: Counterclockwise rotation in radians
    :return: Rotation matrix
    """
    axis = np.asarray(axis)
    axis = axis / math.sqrt(np.dot(axis, axis))
    a = math.cos(theta / 2.0)
    b, c, d = -axis * math.sin(theta / 2.0)
    aa, bb, cc, dd = a * a, b * b, c * c, d * d
    bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
    return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
                     [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
                     [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])


luroth's avatar
.    
luroth committed
117
def create_nband_GeoTiff(file_path, arrays, data_type, keys=None):
luroth's avatar
luroth committed
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
    """Writes geoTiff to filesystem using as many bands as arrays are provided. Optionally annotates metadata per band

    :param file_path: Path to save file
    :param arrays: List ob band arrays
    :param data_type: gdal datatype for all bands
    :param keys: List of dictionaries with metadata per band
    :return:
    """
    driver = gdal.GetDriverByName('GTiff')
    file = file_path

    ds = driver.Create(file, arrays[0].shape[1], arrays[0].shape[0], len(arrays), data_type )

    for i, array in enumerate(arrays):
        band_ = ds.GetRasterBand(i+1)
        band_.WriteArray(array)
luroth's avatar
.    
luroth committed
134
135
136
        if keys is not None:
            for key, value in keys[i].items():
                band_.SetMetadataItem(key, str(value))
luroth's avatar
luroth committed
137
138
139

    ds = None

luroth's avatar
luroth committed
140
141
    return file_path

142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
def read_nband_GeoTiff(file_path):
    """Reads geoTiff from filesystem with as many bands as provided. Optionally reads metadata per band if existing

    :param file_path: Path to save file
    :return: numpy array with bands, annotations
    """
    bands = []
    keys = []

    ds = gdal.Open(file_path)

    for i in range(ds.RasterCount):
        band_ = ds.GetRasterBand(i+1)
        bands.append(band_.ReadAsArray())

        metadata = band_.GetMetadata()
        keys.append(metadata)

    bands = np.stack(bands, axis=2)

    ds = None

    return bands, keys

luroth's avatar
luroth committed
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
def read_nband_mask_GeoTiff(file_path):
    """Reads geoTiff from filesystem with as many bands as provided. Optionally reads metadata per band if existing

        :param file_path: Path to save file
        :return: numpy array with bands, annotations, mask
        """
    bands = []
    keys = []

    ds = gdal.Open(file_path)

    for i in range(ds.RasterCount):
        band_ = ds.GetRasterBand(i + 1)
        bands.append(band_.ReadAsArray())

        metadata = band_.GetMetadata()
        keys.append(metadata)

    bands = np.stack(bands, axis=2)

    mask = band_.GetMaskBand().ReadAsArray()

    ds = None

    return bands, keys, mask

192
193
def create_elevation_diff_GeoTiff(file_path_DTM, file_path_DSM, path_output_file):

luroth's avatar
luroth committed
194
195
    print("Processing", file_path_DSM)

196
    # Read DTM
luroth's avatar
luroth committed
197
198
    dtm_ds = gdal.Open(file_path_DTM)
    dtm_band1 = dtm_ds.GetRasterBand(1)
luroth's avatar
luroth committed
199
    dtm_nodata_val = dtm_band1.GetNoDataValue()
luroth's avatar
luroth committed
200
    dtm_arr = dtm_band1.ReadAsArray()
201
    dtm_transform = dtm_ds.GetGeoTransform()
luroth's avatar
luroth committed
202

203
    # set nondata values to 0
luroth's avatar
luroth committed
204
205
206
    nondata_val = dtm_band1.GetNoDataValue()
    dtm_arr[np.where(dtm_arr == nondata_val)] = 0

207
    # Read DSM
luroth's avatar
luroth committed
208
209
210
    dsm_ds = gdal.Open(file_path_DSM)
    dsm_band1 = dsm_ds.GetRasterBand(1)
    dsm_arr = dsm_band1.ReadAsArray()
211
212
    dsm_transform = dsm_ds.GetGeoTransform()
    # set nondata values to 0
luroth's avatar
luroth committed
213
214
215
    nondata_val = dsm_band1.GetNoDataValue()
    dsm_arr[np.where(dsm_arr == nondata_val)] = 0

luroth's avatar
luroth committed
216
    # test if geoTIFF transforms compatible (same pixel size and orientation - the extent can differ)
217
218
219
220
221
222
223
224
225
226
227
    assert (math.isclose(dtm_transform[1], dsm_transform[1], abs_tol=1e-5) and math.isclose(dtm_transform[2], dsm_transform[2], abs_tol=1e-5)
            and math.isclose(dtm_transform[4], dsm_transform[4], abs_tol=1e-5) and math.isclose(dtm_transform[5], dsm_transform[5], abs_tol=1e-5) )

    # make arrays same size
    diff_x = round((dtm_transform[0] - dsm_transform[0]) / dtm_transform[1])
    diff_x_cut = diff_x if diff_x > 0 else 0
    diff_x_ins = abs(diff_x) if diff_x < 0 else 0
    diff_y = -1 * round((dtm_transform[3] - dsm_transform[3]) / dtm_transform[5])
    diff_y_cut = diff_y if diff_y > 0 else 0
    diff_y_ins = abs(diff_y) if diff_y < 0 else 0

luroth's avatar
luroth committed
228
229
    delta_x = dtm_arr.shape[0] - dsm_arr.shape[0] + diff_x_cut
    delta_x = delta_x if delta_x > 0 else 0
luroth's avatar
luroth committed
230
    delta_y = dtm_arr.shape[1] - dsm_arr.shape[1] + diff_y_cut
luroth's avatar
luroth committed
231
    delta_y = delta_y if delta_y > 0 else 0
232
233


luroth's avatar
luroth committed
234
235
236
237
    dsm_arr_ = dsm_arr

    dsm_arr_ = np.pad(dsm_arr_, ((diff_x_ins, delta_x), (diff_y_ins, delta_y)), 'constant', constant_values= 0)
    dsm_arr_ = dsm_arr_[diff_x_cut:dtm_arr.shape[0] + diff_x_cut, diff_y_cut:dtm_arr.shape[1] + diff_y_cut]
238

luroth's avatar
luroth committed
239
    # difference
240
    ph_arr = np.subtract(dsm_arr_, dtm_arr)
luroth's avatar
luroth committed
241

242
    # convert to float16
luroth's avatar
luroth committed
243
    ph_arr = np.array(ph_arr, dtype=np.float32)
luroth's avatar
luroth committed
244
245
    # set to nadata if no value in height model
    ph_arr[np.where(dsm_arr == nondata_val)] = nondata_val
luroth's avatar
luroth committed
246
    ph_arr[np.where(dtm_arr == dtm_nodata_val)] = nondata_val
luroth's avatar
luroth committed
247
248
249

    # write geoTIFF
    driver = gdal.GetDriverByName('GTiff')
250
    ds_ph = driver.Create(path_output_file, xsize=ph_arr.shape[1], ysize=ph_arr.shape[0], bands=1, eType=gdal.GDT_Float32)
luroth's avatar
luroth committed
251

252
253
    # set geotransform from DTM
    ds_ph.SetGeoTransform(dtm_transform)
luroth's avatar
luroth committed
254
255
256
    # set projection
    ds_ph.SetProjection(dtm_ds.GetProjection())

luroth's avatar
luroth committed
257
258
259
260
261
262
    band1 = ds_ph.GetRasterBand(1)
    band1.WriteArray(ph_arr)
    band1.SetNoDataValue(nondata_val)

    ds_ph.FlushCache()

luroth's avatar
luroth committed
263
264
265
266
    # write to file
    ds_ph = None


luroth's avatar
luroth committed
267
268
269
270
271
272
273
def read_elevation_from_DTM(file_path, points):
    """Reads the elevation of a list of points out of a DTM

    :param file_path: Path to DTM file
    :param points: List of points with x/y coordinates (same projection as the DTM)
    :return: List of elevations
    """
luroth's avatar
luroth committed
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296

    dataset = gdal.Open(file_path)
    band = dataset.GetRasterBand(1)

    cols = dataset.RasterXSize
    rows = dataset.RasterYSize

    transform = dataset.GetGeoTransform()

    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = -transform[5]

    data = band.ReadAsArray(0, 0, cols, rows)

    elevations = []
    for point in points:
        col = int((point[0] - xOrigin) / pixelWidth)
        row = int((yOrigin - point[1]) / pixelHeight)

        elevations.append(data[row][col])

luroth's avatar
.    
luroth committed
297
298
    return elevations

luroth's avatar
luroth committed
299
def zonal_stat(output_path, band_index, raster_path=None, raster_ds=None, shape_str=None, shape_path=None):
luroth's avatar
luroth committed
300
    # Open raster data
luroth's avatar
luroth committed
301
302
303
304
305
    if raster_ds is None:
        raster = gdal.Open(raster_path)
    else:
        raster = raster_ds

luroth's avatar
luroth committed
306
307
308
309
310
311
312
313
314
315
316
317
318
    raster_band = raster.GetRasterBand(band_index)

    # Get raster georeference info
    transform = raster.GetGeoTransform()
    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = transform[5]

    # Container to collect data
    data = []

    # Open geojson polygon file
luroth's avatar
luroth committed
319
320
321
322
323
324
325
    if shape_str:
        samples_polygons_ = geojson.loads(shape_str)
    else:
        with open(shape_path, 'r') as f:
            samples_polygons_ = geojson.load(f)

    samples_polygons = samples_polygons_['features']
luroth's avatar
luroth committed
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353

    # Iterate over polygons, sample each polygon
    for polygon in samples_polygons:
        # Get geometry
        coords = polygon['geometry']['coordinates'][0]
        plot_label = polygon['properties']['plot_label']

        print("Sample no", plot_label)

        # Transform coordinates to path
        sample_path = path.Path(coords)

        # Get extent
        xmin, ymin, xmax, ymax = sample_path.get_extents().extents

        # Specify offset and rows and columns to read
        # Border pixel to ensure plot is in extent
        border = 10
        # Offsets and size
        xoff = math.floor((xmin - xOrigin) / pixelWidth) - border
        yoff = math.floor((yOrigin - ymax) / pixelWidth) - border
        xcount = int((xmax - xmin) / pixelWidth) + 2 * border
        ycount = int((ymax - ymin) / pixelWidth) + 2 * border
        # Recalculate origin
        xmin_inc_border = xoff * pixelWidth + xOrigin
        ymax_inc_border = yOrigin - (yoff * -1 * pixelHeight)

        # Read raster as arrays
luroth's avatar
luroth committed
354
355
        dataraster = None
        i = 1000
luroth's avatar
luroth committed
356
        if (dataraster is None):
luroth's avatar
luroth committed
357
358
            dataraster = raster_band.ReadAsArray(xoff, yoff, xcount, ycount)
            i -= 1
luroth's avatar
luroth committed
359

luroth's avatar
luroth committed
360
361
362
363
        if (dataraster is None):
            print("skipping plot, not readable:",  plot_label)
            continue

luroth's avatar
luroth committed
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
        # Create mask of plot
        # Create coordinate matrix to check if image pixel in plot polygon
        east = np.linspace(xmin_inc_border, xmin_inc_border + dataraster.shape[0] * pixelWidth, dataraster.shape[1])
        north = np.linspace(ymax_inc_border, ymax_inc_border + dataraster.shape[1] * pixelHeight, dataraster.shape[0])
        coords = np.transpose([np.repeat(east, len(north)), np.tile(north, len(east))])

        # Create mask
        sample_mask_image = sample_path.contains_points(coords, radius=abs(pixelWidth))
        sample_mask_image = np.swapaxes(sample_mask_image.reshape(dataraster.shape[1], dataraster.shape[0]), 0, 1)

        # Mask zone of raster
        zoneraster = np.ma.masked_array(dataraster, np.logical_not(sample_mask_image))

        # Calculate statistics of zones
        count = zoneraster.count()
        mean = np.mean(zoneraster)
        std = np.std(zoneraster)
        var = np.var(zoneraster)
        values_percentiles = np.percentile(zoneraster, np.arange(0, 101))

        # add to data container
        value_stat = {'plot_label': plot_label, 'count': count,
                      'mean': mean, 'std': std, 'var': var}
luroth's avatar
luroth committed
387
        values_percentiles = {('percentile_' + f'{i:03d}'): x for i, x in enumerate(values_percentiles)}
luroth's avatar
luroth committed
388
389
390
391
392

        data.append({**value_stat, **values_percentiles})

    data_pd = pd.DataFrame(data, columns=data[0].keys())
    data_pd.to_csv(output_path, index=False)
luroth's avatar
luroth committed
393

luroth's avatar
luroth committed
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
def voronoi_finite_polygons_2d(vor, radius=None):
    """
    Source: http://stackoverflow.com/a/20678647/1595060:

    Reconstruct infinite voronoi regions in a 2D diagram to finite
    regions.

    Parameters
    ----------
    vor : Voronoi
        Input diagram
    radius : float, optional
        Distance to 'points at infinity'.

    Returns
    -------
    regions : list of tuples
        Indices of vertices in each revised Voronoi regions.
    vertices : list of tuples
        Coordinates for revised Voronoi vertices. Same as coordinates
        of input vertices, with 'points at infinity' appended to the
        end.

    """

    if vor.points.shape[1] != 2:
        raise ValueError("Requires 2D input")

    new_regions = []
    new_vertices = vor.vertices.tolist()

    center = vor.points.mean(axis=0)
    if radius is None:
        radius = vor.points.ptp().max()

    # Construct a map containing all ridges for a given point
    all_ridges = {}
    for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
        all_ridges.setdefault(p1, []).append((p2, v1, v2))
        all_ridges.setdefault(p2, []).append((p1, v1, v2))

    # Reconstruct infinite regions
    for p1, region in enumerate(vor.point_region):
        vertices = vor.regions[region]

        if all(v >= 0 for v in vertices):
            # finite region
            new_regions.append(vertices)
            continue

        # reconstruct a non-finite region
        ridges = all_ridges[p1]
        new_region = [v for v in vertices if v >= 0]

        for p2, v1, v2 in ridges:
            if v2 < 0:
                v1, v2 = v2, v1
            if v1 >= 0:
                # finite ridge: already in the region
                continue

            # Compute the missing endpoint of an infinite ridge

            t = vor.points[p2] - vor.points[p1]  # tangent
            t /= np.linalg.norm(t)
            n = np.array([-t[1], t[0]])  # normal

            midpoint = vor.points[[p1, p2]].mean(axis=0)
            direction = np.sign(np.dot(midpoint - center, n)) * n
            far_point = vor.vertices[v2] + direction * radius

            new_region.append(len(new_vertices))
            new_vertices.append(far_point.tolist())

        # sort region counterclockwise
        vs = np.asarray([new_vertices[v] for v in new_region])
        c = vs.mean(axis=0)
        angles = np.arctan2(vs[:, 1] - c[1], vs[:, 0] - c[0])
        new_region = np.array(new_region)[np.argsort(angles)]

        # finish
        new_regions.append(new_region.tolist())

    return new_regions, np.asarray(new_vertices)

def create_voronoi_plot_geoJSON(point_coordinates, plot_size_x, plot_size_y, output_file_path):
    # Convert normalized coordinates to image coordinates
    plant_image_coordinates_x = point_coordinates['x'] * plot_size_x
    plant_image_coordinates_y = point_coordinates['y'] * plot_size_y
    plant_image_coordinates = pd.DataFrame(
        {'x': plant_image_coordinates_x, 'y': plant_image_coordinates_y}, columns=['x', 'y'])

    # Calculate voronoi regions
    vor = Voronoi(plant_image_coordinates)

    # Extract voronoi regions as polygons
    regions, vertices = voronoi_finite_polygons_2d(vor)

    # Clip to plot shape
    box = Polygon([[0, 0], [0, plot_size_y], [plot_size_x, plot_size_y], [plot_size_x, 0]])

    # Create geopandas dataframe for later export as geoJSON
    polygons = [Polygon(vertices[region]).intersection(box) for region in regions]
    gpd_polygons = gpd.GeoDataFrame(geometry=polygons)
    gpd_polygons['id'] = gpd_polygons.index
    gpd_polygons['aera'] = gpd_polygons.area

    # Hirarchical clustering to build groups
    n_points = len(plant_image_coordinates)
    n_points_in_cluser = np.array(range(2, n_points))
    cluster_ns = np.round(n_points / n_points_in_cluser)
    possible_clusters = np.array(np.stack([n_points_in_cluser, cluster_ns], axis=1), dtype=np.int)
    possible_clusters = possible_clusters[possible_clusters[:, 1] > 1]

    for n_points, cluster_n in possible_clusters:
        cluster = AgglomerativeClustering(n_clusters=cluster_n, affinity='euclidean', linkage='ward')
        cluster.fit_predict(plant_image_coordinates[['x', 'y']])
        gpd_polygons['hclust_n' + str(n_points)] = cluster.labels_

    gpd_polygons.to_file(filename=output_file_path, driver='GeoJSON')

luroth's avatar
luroth committed
515

luroth's avatar
.    
luroth committed
516
517
if __name__ == "__main__":
    print("Nothing to run")