process_ph_14m.py 2.67 KB
Newer Older
luroth's avatar
luroth committed
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
from Common import GISfunctions

import os
import shutil
import numpy as np
from scipy import  ndimage
import gdal
import cv2
from pathlib import Path


def process_df(path_dtm, path_dsm, path_output):
    GISfunctions.create_elevation_diff_GeoTiff(path_dtm, path_dsm, path_output)


if __name__ == "__main__":

    path_dtm = r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20181109\14m_M600P\orthoJPEG_DEM.tif'
    path_dtm_smoothed = r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20181109\14m_M600P\orthoJPEG_DEM_mv5median.tif'

    if not Path(path_dtm_smoothed).exists():

        # Read DTM
        print("Reading DTM")
        dtm_ds = gdal.Open(path_dtm)
        dtm_band1 = dtm_ds.GetRasterBand(1)
        dtm_arr = dtm_band1.ReadAsArray()
        dtm_transform = dtm_ds.GetGeoTransform()

        print("Smoothing DTM")
        dtm_smoothed = cv2.medianBlur(dtm_arr, 5)

        print("Saving smoothed DTM")
        driver = gdal.GetDriverByName('GTiff')
        ds_smoothed = driver.Create(path_dtm_smoothed, xsize=dtm_arr.shape[1], ysize=dtm_arr.shape[0], bands=1,
                              eType=gdal.GDT_Float32)
        band1 = ds_smoothed.GetRasterBand(1)
        band1.WriteArray(dtm_smoothed)
        ds_smoothed.SetGeoTransform(dtm_transform)
        ds_smoothed.SetProjection(dtm_ds.GetProjection())
        ds_smoothed = None

    path_dtm = path_dtm_smoothed

    paths_done = [
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190215\14m_M600P',
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190222\14m_M600P',
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190322\14m_M600P'
    ]

    paths = [
        #r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190502\14m_M600P',
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190601\14m_M600P',
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190605\14m_M600P',
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190607\14m_M600P',
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190612\14m_M600P',
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190614\14m_M600P',
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190617\14m_M600P',
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190624\14m_M600P',
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190628\14m_M600P',
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190703\14m_M600P',
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190705\14m_M600P',
        r'O:\UAV\_Processed_\DSP_strickhof_DSWW003\20190708\14m_M600P',
    ]

    for path in paths:

        print("Process", path)

        path_dsm = os.path.join(path, 'orthoJPEG_DEM.tif')
        path_output = os.path.join(path, 'PH_DEM.tif')

        process_df(path_dtm, path_dsm, path_output)