SDII_fig.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

###############################################################
#       DATA ANALYSIS - VERY WET MONSOON SEASONS IN INDIA     #
#      RESULTS PUBLISHED IN GEOPHYSICAL RESEARCH LETTERS      #
#                    Code by Anja Katzenberger                #
###############################################################  
#Link to publication: 
#https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022GL098856

import netCDF4 as nc
import cartopy.crs as ccrs
import cartopy.feature as cf

#import geopandas as gpd
from datetime import date, timedelta

#import geopandas as gpd
import matplotlib.pyplot as plt
#from mpl_toolkits.basemap import Basemap
from matplotlib import cm
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import os


def getdata(file_name):
    x = nc.Dataset(file_name)
    lon = x.variables['lon'][:]
    lat = x.variables['lat'][:]
    pr_kgms2 = x.variables['simple_daily_intensitiy_index_per_time_period'][:]
    pr = pr_kgms2*100 # transform into %
    
    #time = x.variables['time'][:]
    start = date(1850,1,1) # unit of time: number of days since 1850-01-01
    #dates = [start + timedelta(int(n)) for n in time] # transfers dates into more useful format
    
    return(lat,lon,pr)

#dir = '/Users/anjakatzenberger/Dokumente/01_Uni/Masterarbeit/Extremes/Figures/wetdays/netcdfs'
dir = "/p/tmp/anjaka/Extremes/CMIP6_day/SecondaryData/analysis"


files = os.listdir(dir)

for file in files: 
     # if file.endswith("sdii_sub_div.nc") or file.endswith("sdii_sub_div_ensmean.nc"):   
    if file.endswith("sdii_sub_div_ensmean.nc"):      
            (lat,lon,pr) = getdata(dir + '/' + file)
            model = file.split('_')[2]
            if model == "ssp585":
                model = "ensmean"
            
            fig = plt.figure()
            
            ax = plt.axes(projection=ccrs.PlateCarree())
            ax.add_feature(cf.BORDERS)
            #ax.add_feature(cf.LAND)
            ax.coastlines(resolution = '50m')
            gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0, color='gray', alpha=0.3, linestyle='-')
            plt.imshow(pr[0][::-1], interpolation='none', extent = (67,98,6,36),cmap = cm.RdBu)
            
            gl.xlabels_top = False
            gl.ylabels_bottom = True
            gl.ylabels_right = False
            gl.xlocator = mticker.FixedLocator([70,80,90,100])
            gl.ylocator = mticker.FixedLocator([10,20,30,40])
            gl.xformatter = LONGITUDE_FORMATTER
            gl.yformatter = LATITUDE_FORMATTER
            
            plt.title(model)
            
            cbar = plt.colorbar(extend='both')
            plt.clim(-30,30)
            cbar.set_label('%')
            
            plt.xlabel('latitude')
            plt.ylabel('longitude')
            
            
            plt.savefig('/p/tmp/anjaka/Extremes/CMIP6_day/SecondaryData/figures/sdii_' +model + '.png', dpi = 800,bbox_inches='tight')
            plt.show()