Skip to content
Snippets Groups Projects
Select Git revision
  • main default protected
  • optim
2 results
Blame Permalink
main_prepare_inputs.py 5.07 KiB
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 






# -*- coding: UTF-8 -*-
# Python
"""
11-09-2023
@author: jeremy auclair

Main script to run input scripts to generate input data
"""

if __name__ == '__main__':
    
    import sys, os  # system management
    currentdir = os.path.dirname(os.path.realpath(__file__))
    sys.path.insert(0, os.path.dirname(currentdir))
    
    
    from dask.distributed import Client # to parallelise calculations
    import webbrowser  # to open dask dashboard
    from modspa_pixel.config.config import config  # to import config file
    from modspa_pixel.source.code_toolbox import format_duration  # to print memory requirements
    from modspa_pixel.preprocessing.input_toolbox import prepare_directories, set_eodag_config_file
    from modspa_pixel.preprocessing.download_S2 import download_S2_data, extract_zip_archives
    from modspa_pixel.preprocessing.download_ERA5_weather import request_ER5_weather
    from time import time  # to time code excecution
    
    # Get start time of script
    t0 = time()
    
    # Declare paths
    config_file = currentdir + os.sep + 'config' + os.sep + 'config_modspa.json'
    
    # Open config file and load parameters
    config_params = config(config_file)
    
    # Import functions depending on run mode
    mode = config_params.mode
    if mode == 'pixel':
        from modspa_pixel.preprocessing.calculate_ndvi import calculate_ndvi, interpolate_ndvi
    else:
        from modspa_pixel.preprocessing.calculate_ndvi import calculate_ndvi_parcel
        from modspa_pixel.preprocessing.extract_ndvi import extract_ndvi_stats, filter_raw_ndvi, interpolate_ndvi_parcel
        from modspa_pixel.preprocessing.parcel_to_pixel import convert_dataframe_to_xarray
    
    # Run parameters
    max_cpu = config_params.max_cpu
    run_name = config_params.run_name
    preferred_provider = config_params.preferred_provider
    start_date = config_params.start_date
    end_date = config_params.end_date
    ndvi_overwrite = config_params.ndvi_overwrite
    
    # Path parameters
    download_path = config_params.download_path
    shapefile_path = config_params.shapefile_path
    
    # Prepare directories
    prepare_directories(config_file)
    
    # Updating the eodag config file if needed
    set_eodag_config_file(config_params.path_to_config_file, config_params.download_path, config_params.preferred_provider)
    
    # Generate inputs
    
    #===== NDVI =====#
    
    if preferred_provider == 'copernicus':
        ndvi_path = download_path + os.sep + 'SCIHUB' + os.sep + 'NDVI' + os.sep + run_name
    else:
        ndvi_path = download_path + os.sep + 'THEIA' + os.sep + 'NDVI' + os.sep + run_name
    
    # Download optical images
    csv_download_file = ndvi_path + os.sep + 'download.csv'
    raw_images = download_S2_data(start_date, end_date, preferred_provider, csv_download_file, shapefile_path, mode = mode, cloud_cover_limit = config_params.cloud_cover_limit)
    
    # Extract Zip archives
    csv_extract_file = ndvi_path + os.sep + 'extract.csv'
    extracted_images = extract_zip_archives(download_path, raw_images, preferred_provider, csv_extract_file)
    
    # Calculate and interpolate NDVI
    if mode == 'pixel':
        client = Client()
        client
        if config_params.open_browser:
            webbrowser.open('http://127.0.0.1:8787/status', new=2, autoraise=True)

        # Calculate NDVI
        ndvi_precube = calculate_ndvi(extracted_images, config_file)
        
        # Interpolate NDVI
        ndvi_cube = interpolate_ndvi(ndvi_precube, config_file)        
        
        client.close()
    else:
        # Calculate NDVI
        csv_ndvi_path_file = ndvi_path + os.sep + 'ndvi.csv'
        ndvi_path = calculate_ndvi_parcel(extracted_images, ndvi_path, csv_ndvi_path_file, overwrite = ndvi_overwrite, max_cpu = max_cpu)
    
        # Extract NDVI values on shapefile features
        csv_ndvi_extract_file = ndvi_path + os.sep + 'ndvi_extract.csv'
        raw_ndvi = extract_ndvi_stats(ndvi_path, shapefile_path, csv_ndvi_extract_file, max_cpu = 4)
        
        # Filter NDVI values for each feature
        csv_ndvi_filter_file = ndvi_path + os.sep + 'ndvi_filter.csv'
        filtered_ndvi = filter_raw_ndvi(raw_ndvi, csv_ndvi_filter_file, max_cpu = max_cpu)
        
        # Interpolate NDVI values to a daily time step
        csv_ndvi_interp_file = ndvi_path + os.sep + 'ndvi_interp.csv'
        interpolated_ndvi = interpolate_ndvi_parcel(filtered_ndvi, csv_ndvi_interp_file, start_date, end_date, max_cpu = max_cpu)
        
        # Convert DataFrame to netCDF4 file
        ndvi_cube = ndvi_cube_path = ndvi_path + os.sep + run_name + '_NDVI_cube_' + start_date + '_' + end_date + '.nc'
        convert_dataframe_to_xarray(interpolated_ndvi, ndvi_cube, variables = ['NDVI'], data_types = ['u1'])
    
    
    #===== Weather =====#
    
    # Run weather download and formatting script
    weather = request_ER5_weather(config_file, ndvi_cube, shapefile = shapefile_path, mode = mode, raw_S2_image_ref = ndvi_path + os.sep + 'Aurade_test_grid_reference.tif')
    
    # Print formatted runtime
    print('\nExcecution time: ', end = '')
    format_duration(time() - t0)