Mercurial > repos > ecology > landcover_subindicator

comparison land_cover.py @ 0:d79c143be3a7 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyecology/tools-ecology/tree/master/tools/earth commit 4cbace8c3a71b7a1638cfd44a21f5d4b84d2fd15
author ecology
date Mon, 21 Oct 2024 16:46:50 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:d79c143be3a7
1 # Land Cover TE Algorithm Without GEE
2
3 # How to Execute:
4 # - Load the two images (.TIFF) of the rasters
5 # to be processed in the "/data/land_cover/input" folder
6 # - Setting the two filenames in the specific cell of this script (0)
7 # - Run all cells of the script
8 # - Check the script results in the "/data/land_cover/output" folder
9 # - Land Cover Degradation "/data/land_cover/output/lc_dg.tiff"
10 # - Land Cover Transition "/data/land_cover/output/lc_tr.tiff"
11
12 # Librairies import
13 import argparse
14 import json
15 import os
16 import sys
17
18 import matplotlib.pyplot as plt
19
20 import numpy as np
21
22 import rasterio
23 from rasterio.plot import show
24
25 from te_schemas.land_cover import LCLegendNesting
26 from te_schemas.land_cover import LCTransitionDefinitionDeg
27
28 # Methods to manage Rasters
29
30
31 def remap(raster, problem_numbers, alternative_numbers):
32 n_min, n_max = raster.min(), raster.max()
33 replacer = np.arange(n_min, n_max + 1)
34 mask = (problem_numbers >= n_min) & (problem_numbers <= n_max)
35 replacer[problem_numbers[mask] - n_min] = alternative_numbers[mask]
36 raster_replaced = replacer[raster - n_min]
37 return raster_replaced
38
39
40 def saveRaster(dataset, datasetPath, cols, rows, projection, namedataset=None):
41 if namedataset:
42 rasterSet = rasterio.open(datasetPath,
43 'w',
44 driver='GTiff',
45 height=rows,
46 width=cols,
47 count=1,
48 dtype=np.int8,
49 crs=projection,
50 transform=transform, )
51 rasterSet.write(dataset, 1)
52 rasterSet.close()
53 else:
54 rasterSet = rasterio.open(datasetPath,
55 'w',
56 driver='GTiff',
57 height=rows,
58 width=cols,
59 count=1,
60 dtype=np.int8,
61 crs=projection,
62 transform=transform, )
63 rasterSet.write(dataset, 1)
64 rasterSet.set_band_description(1, namedataset)
65 rasterSet.close()
66
67
68 def plot(ndviImage, cmap):
69 src = rasterio.open(ndviImage, 'r')
70 fig, ax = plt.subplots(1, figsize=(12, 12))
71 show(src, cmap=cmap, ax=ax)
72 ax.set_xlabel('Est')
73 ax.set_ylabel('Nord')
74 plt.show()
75
76
77 def plotContour(ndviImage, cmap):
78 src = rasterio.open(ndviImage, 'r')
79 fig, ax = plt.subplots(1, figsize=(12, 12))
80 show(src, cmap=cmap, contour=True, ax=ax)
81 ax.set_xlabel('Est')
82 ax.set_ylabel('Nord')
83 plt.show()
84
85
86 # Setting inputs
87 command_line_args = sys.argv[1:]
88
89
90 parser = argparse.ArgumentParser(description="landcover inputs and outputs")
91 # Add arguments
92 parser.add_argument("--raster_1", help="raster 1")
93 parser.add_argument("--raster_2", help="raster 2")
94 parser.add_argument("--json", help="json")
95
96 args = parser.parse_args(command_line_args)
97
98 # Parse the command line arguments
99
100 # Import data
101
102 path_raster_yi = args.raster_1
103 path_raster_yf = args.raster_2
104 fjson = args.json
105 # Input Rasters
106
107 # Outputs
108 out_dir = os.path.join(os.getcwd(), "data/land_cover/output/")
109 if not os.path.exists(out_dir):
110 os.makedirs(out_dir)
111
112 # Output Rasters
113 path_lc_tr = os.path.join(out_dir, 'lc_tr.tiff')
114 path_lc_bl = os.path.join(out_dir, 'lc_bl.tiff')
115 path_lc_dg = os.path.join(out_dir, 'lc_dg.tiff')
116 path_lc_tg = os.path.join(out_dir, 'lc_tg.tiff')
117 path_change_yf_yi = os.path.join(out_dir, 'change_yf_yi.tiff')
118 path_lc_baseline_esa = os.path.join(out_dir, 'lc_baseline_esa.tiff')
119 path_lc_target_esa = os.path.join(out_dir, 'lc_target_esa.tiff')
120 path_out_img = os.path.join(out_dir, 'stack.tiff')
121
122 # Contours
123 contours_change_yf_yi = os.path.join(out_dir, '/change_yf_yi0.shp')
124
125
126 # Parsing Inputs
127 # Load static inputs
128 # Transition Matrix, ESA Legend, IPCC Legend
129 # Input Raster and Vector Paths
130 params = json.load(open(fjson))
131
132 crs = params.get("crs")
133
134 trans_matrix_val = params.get("trans_matrix")
135 trans_matrix = LCTransitionDefinitionDeg.Schema().load(trans_matrix_val)
136
137 esa_to_custom_nesting = LCLegendNesting.Schema().load(
138 params.get("legend_nesting_esa_to_custom")
139 )
140
141 ipcc_nesting = LCLegendNesting.Schema().load(
142 params.get("legend_nesting_custom_to_ipcc")
143 )
144
145 class_codes = sorted([c.code for c in esa_to_custom_nesting.parent.key])
146 class_positions = [*range(1, len(class_codes) + 1)]
147
148 # Load dynamic inputs
149 # Baseline ESA Raster
150 raster_yi = rasterio.open(path_raster_yi)
151 lc_baseline_esa = raster_yi.read(1)
152 yi_dict_profile = dict(raster_yi.profile)
153
154 for k in yi_dict_profile:
155 print(k.upper(), yi_dict_profile[k])
156
157 # Target ESA Raster
158 raster_yf = rasterio.open(path_raster_yf)
159 lc_target_esa = raster_yf.read(1)
160 yf_dict_profile = dict(raster_yf.profile)
161
162 for k in yf_dict_profile:
163 print(k.upper(), yf_dict_profile[k])
164
165 cols = raster_yi.width
166 rows = raster_yf.height
167 transform = raster_yi.transform
168 projection = raster_yi.crs
169
170 # Setting up output
171 saveRaster(lc_baseline_esa.astype('int8'),
172 path_lc_baseline_esa,
173 cols,
174 rows,
175 projection,
176 "Land_cover_yi")
177
178 saveRaster(lc_target_esa.astype('int8'),
179 path_lc_target_esa,
180 cols,
181 rows,
182 projection,
183 "Land_cover_yf")
184
185 # Algorithm execution
186 # Transition codes are based on the class code indices
187 # (i.e. their order when sorted by class code) -
188 # not the class codes themselves.
189 # So need to reclass the land cover used for the transition calculations
190 # from the raw class codes to the positional indices of those class codes.
191 # And before doing that, need to reclassified initial and
192 # final layers to the IPCC (or custom) classes.
193
194 # Processing baseline raster
195 bl_remap_1 = remap(lc_baseline_esa,
196 np.asarray(esa_to_custom_nesting.get_list()[0]),
197 np.asarray(esa_to_custom_nesting.get_list()[1]))
198 lc_bl = remap(bl_remap_1, np.asarray(class_codes), np.asarray(class_positions))
199
200 saveRaster(lc_bl.astype('int8'),
201 path_lc_bl,
202 cols,
203 rows,
204 projection,
205 "Land_cover_yi")
206
207 # Processing Target Raster
208 tg_remap_1 = remap(lc_target_esa,
209 np.asarray(esa_to_custom_nesting.get_list()[0]),
210 np.asarray(esa_to_custom_nesting.get_list()[1]))
211 lc_tg = remap(tg_remap_1, np.asarray(class_codes), np.asarray(class_positions))
212
213 saveRaster(lc_tg.astype('int8'),
214 path_lc_tg,
215 cols,
216 rows,
217 projection,
218 "Land_cover_yf")
219
220 # Processing Transition Map
221 # Compute transition map (first digit for baseline land cover,
222 # and second digit for target year land cover)
223 lc_tr = (lc_bl * esa_to_custom_nesting.parent.get_multiplier()) + lc_tg
224
225 # Compute land cover transition
226 # Remap persistence classes so they are sequential.
227 # This makes it easier to assign a clear color ramp in QGIS.
228 lc_tr_pre_remap = remap(lc_tr,
229 np.asarray(trans_matrix.get_persistence_list()[0]),
230 np.asarray(trans_matrix.get_persistence_list()[1]))
231
232 saveRaster(lc_tr_pre_remap.astype('int8'),
233 path_lc_tr,
234 cols,
235 rows,
236 projection,
237 "Land_cover_transitions_yi-yf")
238
239 # Compute land cover degradation
240 # definition of land cover transitions as degradation (-1),
241 # improvement (1), or no relevant change (0)
242 lc_dg = remap(lc_tr,
243 np.asarray(trans_matrix.get_list()[0]),
244 np.asarray(trans_matrix.get_list()[1]))
245
246 saveRaster(lc_dg.astype('int8'),
247 path_lc_dg,
248 cols,
249 rows,
250 projection,
251 "Land_cover_degradation")
252
253 # Compute Mutlibands stack
254 # Land Cover Degradation + Baseline ESA
255 # + Target ESA + Land Cover Transition
256 dg_raster = rasterio.open(path_lc_dg, "r")
257 dg = dg_raster.read(1, masked=True)
258 baseline_esa = rasterio.open(path_lc_baseline_esa, "r").read(1, masked=True)
259 target_esa = rasterio.open(path_lc_target_esa, "r").read(1, masked=True)
260 tr = rasterio.open(path_lc_tr, "r").read(1, masked=True)
261
262 band_list = [dg, lc_baseline_esa, lc_target_esa, tr]
263 out_meta = dg_raster.meta.copy()
264 out_meta.update({"count": 4})
265
266 with rasterio.open(path_out_img, 'w', **out_meta) as dest:
267 for band_nr, src in enumerate(band_list, start=1):
268 dest.write(src, band_nr)