Skip to content
Extraits de code Groupes Projets
script_backtracking4shape.py 7,24 ko
Newer Older
  • Learn to ignore specific revisions
  • #!/usr/bin/env python3
    #  -*- coding: utf-8 -*-
    
    """
    Runs 5 days backtracking simulations of oil particles from different sites in Qatar
    
    Author: Thomas Anselain, Earth and Life Institute, UCLouvain, Belgium
    Last modified: 7 April 2022
    """
    #%% Import packages
    
    import sys
    from datetime import timedelta
    from datetime import datetime
    from opendrift.readers import reader_netCDF_CF_generic
    from opendrift.models.openoil import OpenOil
    import os
    from netCDF4 import Dataset
    from mesh_to_landmask import h5_to_landmask
    import calendar
    import numpy as np
    import geopandas as gpd
    
    basedir = "/export/miro/students/tanselain/OpenOil/"
    
    abbr_to_num = {name: num for num, name in enumerate(calendar.month_abbr) if num}
    
    #%% Beginning of loops
    def OpenOil_run(source, mois, year, out) :
        path = basedir + f'Simulation_backward/Shape_desalination_plants/Point_station_buffer_decoupe_point_buffer_{source}.shp'
        df1 = gpd.read_file(path)
        x = len(df1)
        print(x)
        
        num =abbr_to_num[mois] 
        # Take into account the monthly averaged temperature
        b = {'Jan': 22.1, 'Feb': 20.85, 'Mar': 22.25, 'Apr': 24.85, 'May': 25.55, 'Jun': 31.25, 'Jul': 32.6, 'Aug': 33.3, 'Sep': 32.55, 'Oct': 30.9, 'Nov': 28.05, 'Dec': 24.35}
        temp = b[mois]
        
        # Number of loop in the month
        tstart = calendar.timegm([year_int, num, 1, 0, 0, 0])
        tend = calendar.timegm([year_int+1 if num == 12 else year_int, num%12+1, 1, 0, 0, 0])
        ndays = int((tend-tstart)/86400)
        print(mois, ndays)
        
        if mois == 'Sep' :
            nb =16
        else :
            nb = 1
            
        # Simulation loop
        for date in range(nb, ndays+1) :
            #if source == 'umm_al_houl' and mois == 'Jan' : 
            #    continue 
            #if source == 'umm_al_houl' and mois == 'Feb' : 
            #    continue 
            #if source == 'umm_al_houl' and mois == 'Mar' : 
            #    continue 
            #if source == 'umm_al_houl' and mois == 'Apr' : 
            #    continue
            # Find forcing data
            forcdir = basedir + f'Simulation_backward'
            if year == '2016' :
              wind_file = forcdir + f"/Forcing_data/wind/wind_2015_2016.nc" 
            if year == '2017' or year == '2018' : 
              wind_file = forcdir + f"/Forcing_data/wind/wind_2016_2018.nc"
            if year == '2019' or year == '2020' :
              wind_file = forcdir + f"/Forcing_data/wind/wind_2018_2020.nc"
            wave_file = forcdir + f"/Forcing_data/waves/{year}/waves_{year}{num:02d}.nc"
            mercator_file    = forcdir + f"/Forcing_data/currents+tides/{year}/currents+tides_{year}{num:02d}_extended_2px.nc"
            bathy_file = forcdir + "/Forcing_data/bathymetry_persian_gulf.nc"
            
            def add_std_names(filename):
                # to avoid os error
                os.environ["HDF5_USE_FILE_LOCKING"] = "FALSE"
                with Dataset(filename, "r+") as f:
                    for var in f.variables.values():
                        if not hasattr(var, "standard_name"):
                            # ECMWF wind files
                            if var.name == "u10":
                                setattr(var, "standard_name", "x_wind")
                            elif var.name == "v10":
                                setattr(var, "standard_name", "y_wind")
            
            #def invert_ecmwf_wind(filename):
            #     with Dataset(filename, "r+") as f:
            #         latvar = f.variables["latitude"]
            #         lat = np.asarray(latvar[:])
            #         if lat[1] < lat[0]:
            #             latvar[:] = lat[::-1]
            #             uvar, vvar = f.variables["u10"], f.variables["v10"]
            #             for i in range(uvar.shape[0]):
            #                 uvar[i,:] = np.flipud(np.array(uvar[i]))
            #                 vvar[i,:] = np.flipud(np.array(vvar[i]))
            #                
            #invert_ecmwf_wind(wind_file)
            
            
            # add_std_names(wind_file)
            reader_wind = reader_netCDF_CF_generic.Reader(
                filename=wind_file,
                name="wind",
                proj4='+proj=latlong'
            )
              
            
            reader_waves = reader_netCDF_CF_generic.Reader(
                filename=wave_file,
                name='waves',
                proj4="+proj=latlong"
            )
            
            reader_mercator = reader_netCDF_CF_generic.Reader(
                filename=mercator_file,
                name="mercator",
                proj4='+proj=latlong'
            )
        
            reader_bathy = reader_netCDF_CF_generic.Reader(
                filename= bathy_file,
                name="bathy",
                proj4='+proj=latlong'
            )
    
    
    #%% Backward simulations of 5 days in each days of a given month
        
        
            #if year == '2018' and mois == 'Apr' and source == 'ras_laffan' and date <= 3 :
                #continue 
            #if year == '2017' and mois == 'Oct' and source == 'umm_al_houl' and date == 25 :
                #continue 
            o = OpenOil(loglevel=20, weathering_model='noaa')
            o.add_reader([reader_mercator, reader_wind, reader_waves])
            
            o.set_config('general:coastline_action', 'stranding')
            o.set_config('environment:constant:sea_water_temperature', temp)
            o.set_config('environment:constant:sea_water_salinity', 40)
            
            # Seeding some oil particles on 
            oiltype= 'QATAR MARINE' #'ARABIAN HEAVY'  'QATAR MARINE, OIL & GAS' 'QATAR MARINE, PHILLIPS' 'QATAR MARINE'
            start_time=datetime(year_int,num,date)  #2020 ou 2019
            time_set = timedelta(hours=24)
            if mois == 'Dec' and date == 31 :
                time_set = timedelta(hours=23)
            o.seed_from_shapefile( path, number=5*x, layername=None,  #de base 54 et pas 16
                                    featurenum= None, time=start_time+time_set, oil_type=oiltype) #Time begin the next day at midnight
            
            # Adjusting some configuration
            o.set_config('processes:dispersion', True)
            o.set_config('processes:evaporation', True)
            o.set_config('processes:emulsification', True)
            o.set_config('drift:vertical_mixing', True)
            o.set_config('vertical_mixing:timestep', 2)
    
            # Running model
            o.run(steps=5*24*4, time_step=-900, time_step_output=900, outfile=out + f"out_{mois2}_{date}.nc")
                         
            # Print and plot results (additionnal)
            #o.plot_oil_budget()
            #o.plot(fast=False, filename= out + f"out_{mois2}_{date}.png", show=False)
            #o.animation(fast=True) 
            
    
    #%% Give combinations station/month/year assessed
    
    year = sys.argv[1] #(2016, 2017, 2018, 2019, 2020)
    year_int = int(year)
    
    #list_source = ['ras_laffan', 'abu_fontas', 'umm_al_houl', 'ras_laffan_port']
    list_source = ['ras_laffan_port']
    
    list_mois = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']
    #list_mois = ['Sep','Oct','Nov','Dec']
    
    
    for source in list_source :
        for mois in list_mois :
            a = {'Jan': 'january', 'Feb': 'february', 'Mar': 'march', 'Apr': 'april', 'May': 'may', 'Jun': 'june', 'Jul': 'july', 'Aug': 'august', 'Sep': 'september', 'Oct': 'october', 'Nov': 'november', 'Dec': 'december'}
            mois2 = a[mois]
            print(f'The month is {abbr_to_num[mois]:02d}')
            OpenOil_run(source, mois, year, basedir + f"Output_backward/{year}/{source}/{mois2}/")