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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
# -*- coding: UTF-8 -*-
# Python
"""
11-07-2023
@author: jeremy auclair
Adapted : vincent Rivalland
Convert input forcing from xls reference file to spatialized and formated files adapted to run modspa-pixel
"""
import xarray as xr
import os
import pandas as pd
import numpy as np
data_path = './DEV_inputs_test'
size = 10
# Original sets
ndvi_path = data_path + os.sep + 'ndvi_' + str(size) + '.nc'
rain_path = data_path + os.sep + 'rain_' + str(size) + '.tif'
if os.path.isfile(rain_path):
os.remove(rain_path)
ET0_path = data_path + os.sep + 'ET0_' + str(size) + '.tif'
if os.path.isfile(ET0_path):
os.remove(ET0_path)
# Validation sets
xls_NDVI_path = data_path + os.sep + 'xls_NDVI_' + str(size) + '.nc'
xls_Rain_path = data_path + os.sep + 'xls_Rain_' + str(size) + '.tif'
xls_ET0_path = data_path + os.sep + 'xls_ET0_' + str(size) + '.tif'
xls_outputs = data_path + os.sep + 'xls_outputs_' + str(size) + '.nc'
# Modspa excel file - Reference Simmonneaux !
xls_file_path = './SAMIR_xls/SAMIRpixel_Reference_Simonneaux2012.xls'
# %% Lecture de l'ensemble des variables au format dataframe
df_samir_xls = pd.read_excel(xls_file_path, sheet_name=0, header=10, index_col=0)
df_samir_xls = df_samir_xls.loc[:, ~df_samir_xls.columns.str.contains('^Unnamed')]
# df_samir_xls
dates = df_samir_xls.index
# %% Genere NDVI 3D
# Open empty dataset to get structure and reindex with correct dates
empty_dataset = xr.open_dataset(ndvi_path)
empty_dataset = empty_dataset.reindex(time=dates)
# Transpose dimensions
empty_dataset = empty_dataset.transpose('time', 'y', 'x')
# Get the numpy array for 'ndvi'
zero_values = empty_dataset['ndvi'].values
# Transpose the numpy array for 'ndvi'
zero_values = zero_values.transpose([0, 2, 1])
empty_dataset['ndvi'] = empty_dataset.ndvi.transpose('time', 'y', 'x')
# Assign the transposed numpy array back to 'ndvi'
empty_dataset.ndvi.values = zero_values
# Drop ndvi to get empty dataset
empty_dataset = empty_dataset.drop_vars('ndvi')
# Datasets
ndvi_val = empty_dataset.copy(deep=True)
rain_val = empty_dataset.copy(deep=True)
ET0_val = empty_dataset.copy(deep=True)
outputs_val = empty_dataset.copy(deep=True)
# Inputs
ndvi_val['NDVI'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype='uint8'))
rain_val['Rain'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype='uint16'))
ET0_val['ET0'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype='uint16'))
# Outputs
liste_variables = df_samir_xls.columns[0:-2]
# for var in val_results.columns:
for var in liste_variables:
outputs_val[var] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype='int16'))
for x in ndvi_val.coords['x'].values:
for y in ndvi_val.coords['y'].values:
# Inputs
ndvi_val['NDVI'].loc[{'x': x, 'y': y}] = np.round(df_samir_xls['NDVI'].values * 255)
rain_val['Rain'].loc[{'x': x, 'y': y}] = np.round(df_samir_xls['Rain'].values * 1000)
ET0_val['ET0'].loc[{'x': x, 'y': y}] = np.round(df_samir_xls['ET0'].values * 1000)
# Outputs
# for var in list(outputs_val.keys()):
for var in liste_variables:
outputs_val[var].loc[{'x': x, 'y': y}] = np.round(df_samir_xls[var].values * 100)
# Add precip
outputs_val['Rain'] = rain_val['Rain'].copy(deep=True) / 10
# Reorganize dimension order
ndvi_val = ndvi_val.transpose('time', 'y', 'x')
rain_val = rain_val.transpose('time', 'y', 'x')
ET0_val = ET0_val.transpose('time', 'y', 'x')
# Save datasets
# Inputs
ndvi_val.to_netcdf(xls_NDVI_path, encoding={"NDVI": {"dtype": "u1", "_FillValue": 0}})
rain_val.Rain.rio.to_raster(xls_Rain_path, dtype='uint16')
ET0_val.ET0.rio.to_raster(xls_ET0_path, dtype='uint16')
# Create encoding dictionnary
for variable in list(outputs_val.keys()):
# Write encoding dict
encoding_dict = {}
encod = {}
encod['dtype'] = 'i2'
encoding_dict[variable] = encod
# Outputs
outputs_val.to_netcdf(xls_outputs, encoding=encoding_dict)