Skip to content

Latest commit

 

History

History
140 lines (114 loc) · 3.58 KB

test_drought.md

File metadata and controls

140 lines (114 loc) · 3.58 KB
jupyter
jupytext kernelspec
formats text_representation
ipynb,md
extension format_name format_version jupytext_version
.md
markdown
1.3
1.16.6
display_name language name
climate_env
python
climate_env
 
import os
import sys
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import argparse
import geojson

# import drought indicies
from climate_drought import drought_indices as dri, config

INDEX_MAP = {
    'SPI_ECMWF': dri.SPI_ECMWF,
    'SPI_GDO': dri.SPI_GDO,
    'SPI_NCG': dri.SPI_NCG,
    'SMA_ECMWF': dri.SMA_ECMWF,
    'SMA_GDO': dri.SMA_GDO,
    'fAPAR': dri.FPAR_GDO,
    'CDI': dri.CDI,
    'FEATURE_SAFE': dri.FEATURE_SAFE,
    'UTCI': dri.UTCI
}
class Drought:
    def __init__(self,product,latitude,longitude,start_date,end_date):
        # Setup paramaters
        self.verbose = True

        self.indir = './input'
        self.outdir = './output'
        self.oformat = 'GeoJSON'
        self.product = product
        self.start_date = start_date
        self.end_date = end_date
        
        # Convert latitude and longitude strings to lists
        self.latitude = [float(item) for item in latitude.replace('[','').replace(']','').split(',')]
        self.longitude = [float(item) for item in longitude.replace('[','').replace(']','').split(',')]

        # setup config and args
        self.cfg = config.Config(self.outdir,self.indir)
        self.args = config.AnalysisArgs(latitude,longitude,
            start_date,end_date,product=product,oformat=self.oformat)
product = "UTCI"
latitude = '52.5' 
longitude = '1.25'
start_date = '20220101'
end_date = '20241231'

print("Running {} for {} {} from {} to {}".format(product, 
    latitude, longitude, start_date, end_date))
 
obj = Drought(product,latitude,longitude,start_date,end_date)
# Run processing
def drought_index(obj) -> dri.DroughtIndex:
   return INDEX_MAP[obj.product](obj.cfg, obj.args)

idx = drought_index(obj)

print("Computing {} index for {} to {}.".format(obj.product, 
    obj.cfg.baseline_start, obj.cfg.baseline_end))

if os.path.exists(idx.output_file_path):
    print("Processed file '{}' already exists.".format(idx.output_file_path))
else:
    idx.download()
    idx.process()
    print("Downloading and processing complete for '{}' completed with format {}.".format(idx.output_file_path, obj.oformat))

if os.path.exists(idx.output_file_path):
    exit_code = 1
    print("{} processing complete, generated {}".format(product, idx.output_file_path))

else:
    print("Processing failed, {} does not exist".format(idx.output_file_path))
import geopandas as gpd

# Load in data and display then plot
df = gpd.read_file(idx.output_file_path)
print(df)
# plotting

fig, ax1 = plt.subplots()
ax1.plot(df._date,df.spi,color='b',label='spi')
ax1.set_ylabel('SPI [blue]')
tick_list = df._date.values[::3]
plt.xticks(rotation=45, ticks=tick_list)
minv = np.nanmin([np.nanmin(df.hindex),np.nanmin(df.spi)])
maxv = np.nanmax([np.nanmax(df.hindex),np.nanmax(df.spi)])
if product == 'UTCI':
    ax1.plot(df._date,df.hindex,color='g',label='utci')
    ax1.set_ylabel('SPI [blue], Health index [green]')
    ax2 = ax1.twinx()
    ax2.plot(df._date, df.utci, color = 'r', label = 'utci')
    ax2.set_ylabel('UTCI [degC, red]')
    ax1.fill_between(df._date,1.5,maxv, 
                     color = 'C1', alpha=0.3, interpolate = True)
    ax1.fill_between(df._date,-1.5,minv, 
                     color = 'C0', alpha=0.3, interpolate = True)
plt.tight_layout()
plt.show()