Skip to content

Commit 1278973

Browse files
sdunesmesdunesme
committed
Update workflow
1 parent 2f4ddf7 commit 1278973

5 files changed

Lines changed: 306 additions & 1109 deletions

File tree

‎requirements/fct.txt‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,4 +14,5 @@ duecredit==0.9.2
1414
python-dotenv==0.21.1
1515
matplotlib==3.6.3
1616
rastachimp==0.3
17-
shapelysmooth==0.2.0
17+
shapelysmooth==0.2.0
18+
opencv-python==4.10.0.84

‎tutorials/01-drainage_plan/01-drainage.py‎

Lines changed: 104 additions & 80 deletions
Original file line numberDiff line numberDiff line change
@@ -1,128 +1,144 @@
1-
# coding: utf-8
2-
3-
'''
4-
Goals of step 01-drainage :
5-
- Get a clean Flow Accumulation raster in km2
6-
-
7-
'''
1+
import os
2+
from fct.config import config
3+
from multiprocessing import cpu_count
4+
p = cpu_count()/2
85

96
# Create your tileset
10-
117
from fct.cli import Tiles
12-
Tiles.CreateTileset('dem', 10000.0,
13-
tileset1 = '/data/sdunesme/fct/tests_1m/fct_workdir/10k_tileset.gpkg',
14-
tileset2 = '/data/sdunesme/fct/tests_1m/fct_workdir/10kbis_tileset.gpkg')
8+
Tiles.CreateTileset('dem5m', 30000.0,
9+
tileset1 = os.path.join(config.workdir, 'default_tileset.gpkg'),
10+
tileset2 = os.path.join(config.workdir, 'aux_tileset.gpkg'))
1511

16-
# Prepare the DEM tiles and VRT
12+
# ATTENTION : Bien vérifier que les tilesets couvrent bien tout le DEM avant de continuer
13+
14+
# Copy the Hydrographic Reference to outputs/GLOBAL/REFHYDRO
15+
if not os.path.isdir(os.path.join(config.workdir, 'GLOBAL/REFHYDRO/')):
16+
os.makedirs(os.path.join(config.workdir, 'GLOBAL/REFHYDRO/'))
1717

18+
# COPY input REFHYDRO to /local/sdunesme/run_rmc_5m/GLOBAL/REFHYDRO/REFERENTIEL_HYDRO.shp
19+
20+
# Prepare the DEM tiles and VRT
1821
from fct.cli import Tiles
19-
Tiles.DatasourceToTiles('dem', '10k', 'dem', processes=32)
20-
Tiles.DatasourceToTiles('dem', '10kbis', 'dem', processes=32)
22+
Tiles.DatasourceToTiles('dem5m', 'default', 'dem', processes=p) # Environ 120G pour 32*20km
23+
Tiles.DatasourceToTiles('dem5m', 'aux', 'dem', processes=p)
2124

2225
from fct.tileio import buildvrt
23-
buildvrt('10k', 'dem')
24-
buildvrt('10kbis', 'dem')
26+
buildvrt('default', 'dem')
27+
buildvrt('aux', 'dem')
2528

26-
# If you have two different scales DEM, you can fill the precise one with the less precise
2729
# First step when you have only one DEM : Smoothing
2830
from fct.drainage import PrepareDEM
2931
params = PrepareDEM.SmoothingParameters()
30-
params.window=15
32+
#params.window=15 # Uncomment for 1m resolution DEM
3133

32-
PrepareDEM.MeanFilter(params, overwrite=True, processes=32, tileset='10k')
33-
PrepareDEM.MeanFilter(params, overwrite=True, processes=32, tileset='10kbis')
34+
PrepareDEM.MeanFilter(params, overwrite=True, processes=p, tileset='default')
35+
PrepareDEM.MeanFilter(params, overwrite=True, processes=p, tileset='aux')
36+
37+
from fct.tileio import buildvrt
38+
buildvrt('default', 'smoothed')
39+
buildvrt('aux', 'smoothed')
40+
41+
# Drape hydrography network
42+
from fct.drainage import Drape
43+
params = Drape.Parameters()
44+
#params.elevations = 'smoothed'
45+
params.stream_network = 'network-cartography-ready'
46+
params.draped = 'stream-network-draped'
47+
Drape.DrapeNetworkAndAdjustElevations(params)
48+
Drape.SplitStreamNetworkIntoTiles(params, tileset='default')
49+
Drape.SplitStreamNetworkIntoTiles(params, tileset='aux')
3450

3551
# Fill sinks
3652
from fct.drainage import DepressionFill
3753
params = DepressionFill.Parameters()
38-
params.elevations = 'smoothed'
39-
params.offset = -1.0 # burn offset in meters
54+
#params.elevations = 'smoothed'
55+
params.offset = 1.0 # burn offset in meters
4056
params.exterior_data = 9000.0 # exterior value
4157

42-
DepressionFill.LabelWatersheds(params, overwrite=True, processes=32)
43-
DepressionFill.LabelWatersheds(params, overwrite=True, processes=32, tileset='10kbis')
58+
DepressionFill.LabelWatersheds(params, overwrite=True, processes=p)
59+
DepressionFill.LabelWatersheds(params, overwrite=True, processes=p, tileset='aux')
4460

4561
from fct.tileio import buildvrt
46-
buildvrt('10k', 'dem-watershed-labels')
47-
buildvrt('10kbis', 'dem-watershed-labels')
62+
buildvrt('default', 'dem-watershed-labels')
63+
buildvrt('aux', 'dem-watershed-labels')
4864

4965
DepressionFill.ResolveWatershedSpillover(params, overwrite=True)
50-
DepressionFill.ResolveWatershedSpillover(params, overwrite=True, tileset='10kbis')
66+
DepressionFill.ResolveWatershedSpillover(params, overwrite=True, tileset='aux')
5167

52-
DepressionFill.DispatchWatershedMinimumZ(params, overwrite=True, processes=32)
53-
DepressionFill.DispatchWatershedMinimumZ(params, overwrite=True, processes=32, tileset='10kbis')
68+
DepressionFill.DispatchWatershedMinimumZ(params, overwrite=True, processes=p)
69+
DepressionFill.DispatchWatershedMinimumZ(params, overwrite=True, processes=p, tileset='aux')
5470

5571
from fct.tileio import buildvrt
56-
buildvrt('10k', 'dem-filled-resolved')
57-
buildvrt('10kbis', 'dem-filled-resolved')
72+
buildvrt('default', 'dem-filled-resolved')
73+
buildvrt('aux', 'dem-filled-resolved')
5874

5975
# Resolve flats
6076
from fct.drainage import BorderFlats
6177
params = BorderFlats.Parameters()
62-
BorderFlats.LabelBorderFlats(params=params, processes=32)
63-
BorderFlats.LabelBorderFlats(params=params, processes=32, tileset='10kbis')
78+
BorderFlats.LabelBorderFlats(params=params, processes=p)
79+
BorderFlats.LabelBorderFlats(params=params, processes=p, tileset='aux')
6480

6581
from fct.tileio import buildvrt
66-
buildvrt('10k', 'dem-flat-labels')
67-
buildvrt('10kbis', 'dem-flat-labels')
82+
buildvrt('default', 'dem-flat-labels')
83+
buildvrt('aux', 'dem-flat-labels')
6884

6985
BorderFlats.ResolveFlatSpillover(params=params)
70-
BorderFlats.ResolveFlatSpillover(params=params, tileset='10kbis')
86+
BorderFlats.ResolveFlatSpillover(params=params, tileset='aux')
7187

72-
BorderFlats.DispatchFlatMinimumZ(params=params, overwrite=True, processes=32)
73-
BorderFlats.DispatchFlatMinimumZ(params=params, overwrite=True, processes=32, tileset='10kbis')
88+
BorderFlats.DispatchFlatMinimumZ(params=params, overwrite=True, processes=p)
89+
BorderFlats.DispatchFlatMinimumZ(params=params, overwrite=True, processes=p, tileset='aux')
7490

7591
from fct.tileio import buildvrt
76-
buildvrt('10k', 'dem-drainage-resolved')
77-
buildvrt('10kbis', 'dem-drainage-resolved')
92+
buildvrt('default', 'dem-drainage-resolved')
93+
buildvrt('aux', 'dem-drainage-resolved')
7894

7995
# FlatMap.DepressionDepthMap is useful if you want to check which flat areas have been resolved
8096

8197
# Flow direction
8298
from fct.drainage import FlowDirection
8399
params = FlowDirection.Parameters()
84100
params.exterior = 'off'
85-
FlowDirection.FlowDirection(params=params, overwrite=True, processes=32)
86-
FlowDirection.FlowDirection(params=params, overwrite=True, processes=32, tileset='10kbis')
101+
FlowDirection.FlowDirection(params=params, overwrite=True, processes=p)
102+
FlowDirection.FlowDirection(params=params, overwrite=True, processes=p, tileset='aux')
87103

88104
from fct.tileio import buildvrt
89-
buildvrt('10k', 'flow')
90-
buildvrt('10kbis', 'flow')
105+
buildvrt('default', 'flow')
106+
buildvrt('aux', 'flow')
91107

92108
# Flow tiles inlets/outlets graph
93109
from fct.drainage import Accumulate
94110
params = Accumulate.Parameters()
95111
params.elevations = 'dem-drainage-resolved'
96112

97-
Accumulate.Outlets(params=params, processes=32)
98-
Accumulate.Outlets(params=params, processes=32, tileset='10kbis')
113+
Accumulate.Outlets(params=params, processes=p)
114+
Accumulate.Outlets(params=params, processes=p, tileset='aux')
99115

100116
Accumulate.AggregateOutlets(params)
101-
Accumulate.AggregateOutlets(params, tileset='10kbis')
117+
Accumulate.AggregateOutlets(params, tileset='aux')
102118

103119
# Resolve inlets/outlets graph
104120
Accumulate.InletAreas(params=params)
105-
Accumulate.InletAreas(params=params, tileset='10kbis')
121+
Accumulate.InletAreas(params=params, tileset='aux')
106122

107123
# Flow accumulation
108-
Accumulate.FlowAccumulation(params=params, overwrite=True, processes=32)
109-
Accumulate.FlowAccumulation(params=params, overwrite=True, processes=32, tileset='10kbis')
124+
Accumulate.FlowAccumulation(params=params, overwrite=True, processes=p)
125+
Accumulate.FlowAccumulation(params=params, overwrite=True, processes=p, tileset='aux')
110126

111127
from fct.tileio import buildvrt
112-
buildvrt('10k', 'acc')
113-
buildvrt('10kbis', 'acc')
128+
buildvrt('default', 'acc')
129+
buildvrt('aux', 'acc')
114130

115131
# Stream Network from sources
116132
from fct.drainage import StreamSources
117133

118134
StreamSources.InletSources(params)
119-
StreamSources.InletSources(params, tileset='10kbis')
135+
StreamSources.InletSources(params, tileset='aux')
120136

121-
StreamSources.StreamToFeatureFromSources(min_drainage=500, processes=32)
122-
StreamSources.StreamToFeatureFromSources(min_drainage=500, processes=32, tileset='10kbis')
137+
StreamSources.StreamToFeatureFromSources(min_drainage=500, processes=p)
138+
StreamSources.StreamToFeatureFromSources(min_drainage=500, processes=p, tileset='aux')
123139

124140
StreamSources.AggregateStreamsFromSources()
125-
StreamSources.AggregateStreamsFromSources(tileset='10kbis')
141+
StreamSources.AggregateStreamsFromSources(tileset='aux')
126142

127143
# Find NoFlow pixels from RHTS
128144
from fct.drainage import FixNoFlow
@@ -131,53 +147,61 @@
131147
params.noflow = 'noflow'
132148
params.fixed = 'noflow-targets'
133149

134-
FixNoFlow.DrainageRaster(params=params, processes=32)
135-
FixNoFlow.DrainageRaster(params=params, processes=32, tileset='10kbis')
150+
FixNoFlow.DrainageRaster(params=params, processes=p)
151+
FixNoFlow.DrainageRaster(params=params, processes=p, tileset='aux')
136152

137153
from fct.tileio import buildvrt
138-
buildvrt('10k', 'drainage-raster-from-sources')
139-
buildvrt('10kbis', 'drainage-raster-from-sources')
154+
buildvrt('default', 'drainage-raster-from-sources')
155+
buildvrt('aux', 'drainage-raster-from-sources')
140156

141-
FixNoFlow.NoFlowPixels(params=params, processes=32)
142-
FixNoFlow.NoFlowPixels(params=params, processes=32, tileset='10kbis')
157+
FixNoFlow.NoFlowPixels(params=params, processes=p)
158+
FixNoFlow.NoFlowPixels(params=params, processes=p, tileset='aux')
143159

144160
FixNoFlow.AggregateNoFlowPixels(params)
145-
FixNoFlow.AggregateNoFlowPixels(params, tileset='10kbis')
161+
FixNoFlow.AggregateNoFlowPixels(params, tileset='aux')
146162

147163
# Fix NoFlow (create TARGETS shapefile and fix Flow Direction data)
148-
FixNoFlow.FixNoFlow(params, tileset1='10k', tileset2='10kbis', fix=True)
164+
FixNoFlow.FixNoFlow(params, tileset1='default', tileset2='aux', fix=True)
165+
FixNoFlow.FixNoFlow(params, tileset1='aux', tileset2='default', fix=True)
149166

150167
# Remake FlowAccumulation with NoFlow pixels fixed
168+
# Il est possible de déplacer un peu les sources à la main et relancer à partir d'ici pour recalculer le RHTS
151169
from fct.drainage import Accumulate
152170
params = Accumulate.Parameters()
153171
params.elevations = 'dem-drainage-resolved'
172+
# params.exterior_flow = 'exterior-inlet' # Uncomment this if you have exterior EXTERIOR INLET points
154173

155-
Accumulate.Outlets(params=params, processes=32)
156-
Accumulate.AggregateOutlets(params)
157-
Accumulate.InletAreas(params=params)
158-
Accumulate.FlowAccumulation(params=params, overwrite=True, processes=32)
174+
Accumulate.Outlets(params=params, processes=p, tileset='default')
175+
Accumulate.AggregateOutlets(params, tileset='default')
176+
Accumulate.InletAreas(params=params, tileset='default')
177+
Accumulate.FlowAccumulation(params=params, overwrite=True, processes=p, tileset='default')
159178

160179
# Remake stream Network from sources with NoFlow pixels fixed
180+
from fct.drainage import Accumulate
181+
params = Accumulate.Parameters()
182+
params.elevations = 'dem-drainage-resolved'
183+
# params.exterior_flow = 'exterior-inlet' # Uncomment this if you have exterior EXTERIOR INLET points
161184
from fct.drainage import StreamSources
162185

163-
StreamSources.InletSources(params)
164-
StreamSources.StreamToFeatureFromSources(min_drainage=500, processes=32)
165-
StreamSources.AggregateStreamsFromSources()
186+
StreamSources.InletSources(params, tileset='default')
187+
StreamSources.StreamToFeatureFromSources(min_drainage=500, processes=p, tileset='default')
188+
StreamSources.AggregateStreamsFromSources(tileset='default')
189+
190+
# IF NETWORK IS DISCONNECTED, TRY TO RESTART THE WORKFLOW FROM "Remake flow accumulation" WITH THE OTHER TILESET
166191

167192
# Identify network nodes
168193
from fct.drainage import IdentifyNetworkNodes
169194
params = IdentifyNetworkNodes.Parameters()
170195

171-
IdentifyNetworkNodes.IdentifyNetworkNodes(params)
172-
IdentifyNetworkNodes.JoinSourcesAttributes(params)
173-
174-
# Sur le réseau hydro => rang de Hack => Créer champ AXIS = 16 derniers caractères de liens_vers_coursd_eau => extraction des sources
196+
IdentifyNetworkNodes.IdentifyNetworkNodes(params, tileset='default')
197+
IdentifyNetworkNodes.JoinSourcesAttributes(params, tileset='default')
175198

199+
# Join network attributes and aggregate by axis
176200
import os
177-
if not os.path.isdir('/data/sdunesme/fct/tests_1m/fct_workdir/GLOBAL/MEASURE'):
178-
os.mkdir('/data/sdunesme/fct/tests_1m/fct_workdir/GLOBAL/MEASURE')
201+
if not os.path.isdir(f'{config.workdir}/GLOBAL/MEASURE'):
202+
os.mkdir(f'{config.workdir}/GLOBAL/MEASURE')
179203

180204
from fct.drainage import JoinNetworkAttributes
181-
JoinNetworkAttributes.JoinNetworkAttributes('/data/sdunesme/fct/tests_1m/fct_workdir/GLOBAL/DEM/SOURCES_IDENTIFIED_10K.shp', '/data/sdunesme/fct/tests_1m/fct_workdir/GLOBAL/DEM/NETWORK_IDENTIFIED_10K.shp', '/data/sdunesme/fct/tests_1m/fct_workdir/GLOBAL/DEM/RHTS.shp')
182-
# JoinNetworkAttributes.UpdateLengthOrder('/data/sdunesme/fct/tests_1m/fct_workdir/GLOBAL/DEM/RHTS.shp', '/data/sdunesme/fct/tests_1m/fct_workdir/GLOBAL/DEM/RHTS.shp')
183-
JoinNetworkAttributes.AggregateByAxis('/data/sdunesme/fct/tests_1m/fct_workdir/GLOBAL/DEM/RHTS.shp', '/data/sdunesme/fct/tests_1m/fct_workdir/GLOBAL/MEASURE/REFAXIS.shp')
205+
JoinNetworkAttributes.JoinNetworkAttributes(f'{config.workdir}/GLOBAL/DEM/SOURCES_IDENTIFIED_DEFAULT.shp', f'{config.workdir}/GLOBAL/DEM/NETWORK_IDENTIFIED_DEFAULT.shp', f'{config.workdir}/GLOBAL/DEM/RHTS.shp')
206+
# JoinNetworkAttributes.UpdateLengthOrder(f'{config.workdir}/GLOBAL/DEM/RHTS.shp', f'{config.workdir}/GLOBAL/DEM/RHTS.shp')
207+
JoinNetworkAttributes.AggregateByAxis(f'{config.workdir}/GLOBAL/DEM/RHTS.shp', f'{config.workdir}/GLOBAL/MEASURE/REFAXIS.shp')

0 commit comments

Comments
 (0)