Commit f12ac32c authored by Jared Rennie's avatar Jared Rennie
Browse files

initial commit

parents
################################################################################
# Written By Jared Rennie
################################################################################
# Import Packages
import sys, time, datetime, os
import numpy as np
import numpy.ma as ma
import geopandas as gpd
import pandas as pd
import xarray as xr
import calendar
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
from scipy import stats
import warnings
warnings.filterwarnings("ignore")
# Define directories
main_directory="/store/sfcnet/jared/cpo_drought"
data_directory=main_directory+'/input_data/eddi'
shapefile_directory=main_directory+'/results_shapefile'
csv_directory=main_directory+'/results_csv'
plot_directory=main_directory+'/results_plots'
#################################################
# BEGIN PROGRAM
start=time.time()
# Read in Arguments
if len(sys.argv) < 3:
sys.exit("USAGE: <Year-Month> <Month_Lag>\n example: ipython plot_eddi.py 201608 01")
inDate=sys.argv[1]
inLag=sys.argv[2]
inCode1='eddi_'+inLag+'-raw'
inCode2='eddi_'+inLag+'-cat'
inCode3='eddi_'+inLag+'-usdm'
if inLag=='01':
outElement='(1-Month)'
if inLag=='03':
outElement='(3-Month)'
if inLag=='06':
outElement='(6-Month)'
if inLag=='12':
outElement='(12-Month)'
lastDay=str(calendar.monthrange(int(inDate[0:4]), int(inDate[4:6]))[1])
fullDate=inDate+lastDay
#################################################
# READ IN SHAPEFILE
input_shapefile=main_directory+'/input_shapefile/cb_2018_us_county_500k.shp'
print("READ IN SHAPEFILE: ",input_shapefile)
geo_shapefile = gpd.read_file(input_shapefile)
geo_shapefile['GEOID2_INT'] = np.array(geo_shapefile['GEOID'].values,dtype='i')
# Get Projection
projection=geo_shapefile.crs
#################################################
# READ IN GRID
num_lons=464
num_lats=224
min_lon=-125
min_lat=25
cell_size=0.125
missing=-9999.
max_lon=-67-cell_size
max_lat=53-cell_size
lons=np.linspace(min_lon,max_lon,num_lons)
lats=np.linspace(min_lat,max_lat,num_lats)
lats=np.flip(lats)
input_grid=data_directory+'/EDDI_ETrs_'+inLag+'mn_'+fullDate+'.asc'
print("READING IN DATA: "+str(input_grid))
file_handle = open(input_grid, 'r')
contents = file_handle.readlines()
file_handle.close()
# create an Empty DataFrame object with fake row
data_month=pd.DataFrame()
data_month['lat'] = [missing]
data_month['lon'] = [missing]
data_month[inCode1] = [missing]
# Run Through Data
lat_counter=0
for line_counter in range(len(contents)):
if line_counter > 5:
#print(line_counter)
latitude=lats[lat_counter]
line=contents[line_counter].strip().split()
lon_counter=0
for val in line:
longitude=lons[lon_counter]
if float(val) > missing:
data_month = data_month.append({'lat' : latitude, 'lon' : longitude, inCode1 : float(val)}, ignore_index = True)
lon_counter+=1
lat_counter+=1
# Remove Fake Row
data_month=data_month.drop(data_month.index[0])
# Convert table from Pandas to GeoPandas
geo_grid=gpd.GeoDataFrame(data_month,geometry=gpd.points_from_xy(data_month.lon, data_month.lat))
geo_grid['geometry']=geo_grid['geometry'].buffer(cell_size/2.)
# Set Projection to same as Shapefile
geo_grid=geo_grid.set_crs(projection)
#################################################
# Spatial Join
print("SPATIAL JOIN")
# Perform the Join
result=gpd.sjoin(geo_shapefile,geo_grid,how="left", op="intersects")
# Take the Aggregate Mean By County
result_mean=result.groupby(['GEOID'],as_index=False).mean()[['GEOID', inCode1]]
result_mean[inCode1]=result_mean[inCode1].round(decimals=2)
result_mean=result_mean.dropna()
# Join Results with CONUS Counties, to reconcile missing counties during Spatial Join
conus_counties=pd.read_csv(main_directory+'/counties.csv')
conus_counties['GEOID']=conus_counties['GEOID'].astype(str).apply(lambda x: x.zfill(5))
merge_conus=conus_counties.merge(result_mean, on='GEOID', how='left')
result_mean=merge_conus
# Insert and Reformat Columns
result_mean[inCode1]=result_mean[inCode1].fillna(value=-9999.)
result_mean.insert(0, "Date", np.full(result_mean.shape[0],inDate), True)
result_mean.insert(3, inCode2, np.full(result_mean.shape[0],'-1'), True)
result_mean.insert(4, inCode3, np.full(result_mean.shape[0],'-1'), True)
# Run through each county and Create USDM drought category.
for index, row in result_mean.iterrows():
val=row[inCode1]
if val==-9999.:
result_mean.at[index,inCode2]='N/A'
result_mean.at[index,inCode3]='N/A'
continue
cat1='-1'
cat2='-1'
if val>=0.5:
cat1='ED0'
cat2='0'
if val>=0.8:
cat1='ED1'
cat2='1'
if val>=1.3:
cat1='ED2'
cat2='2'
if val>=1.6:
cat1='ED3'
cat2='3'
if val>=2.0:
cat1='ED4'
cat2='4'
if val > -0.5 and val < 0.5:
cat1='-1'
cat2='-1'
if val<=-0.5:
cat1='EW0'
cat2='-1'
if val<=-0.8:
cat1='EW1'
cat2='-1'
if val<=-1.3:
cat1='EW2'
cat2='-1'
if val<=-1.6:
cat1='EW3'
cat2='-1'
if val<=-2.0:
cat1='EW4'
cat2='-1'
result_mean.at[index,inCode2]=cat1
result_mean.at[index,inCode3]=cat2
# Create new shapfeile of Results
out_shapefile=geo_shapefile.merge(result_mean, on='GEOID')
# Save as new shapefile
finalShapeFile=shapefile_directory+"/eddi-"+inLag+"_"+inDate+".shp"
print("OUTPUT TO: "+str(finalShapeFile))
out_shapefile.to_file(finalShapeFile)
# Fix SD FIPS so it can be 46102 in shapefile (to plot), but 46113 in csv (to match with health data)
fipsChange=result_mean.index[result_mean['GEOID'] == '46102'].tolist()[0]
result_mean.at[fipsChange,'GEOID']='46113'
# Save as new csv
finalTextFile=csv_directory+"/eddi-"+inLag+"_"+inDate+".csv"
print("OUTPUT TO CSV: "+str(finalTextFile))
result_mean.to_csv(finalTextFile,index=False)
#################################################
# PLOTTING RAW VALUE
print("PLOTTING RAW VALUE")
dpi=300
# Set Bounds
minLat = 22
maxLat = 50
minLon = -120
maxLon = -73
# Set Up Figure
fig= plt.figure(num=1, figsize=(8,5), dpi=dpi, facecolor='w', edgecolor='k')
ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.LambertConformal())
ax.set_extent([minLon, maxLon, minLat, maxLat], crs=ccrs.Geodetic())
plt.style.use('dark_background')
ax.set_facecolor('black')
# Plot Data For Each County
attribute_counter=0
for county in shpreader.Reader(finalShapeFile).geometries():
val=out_shapefile.iloc[attribute_counter][inCode2]
if val != 'N/A':
if val=='EW4':
outColor='#0000ff'
if val=='EW3':
outColor='#4169e1'
if val=='EW2':
outColor='#1d90ff'
if val=='EW1':
outColor='#00bfff'
if val=='EW0':
outColor='#8ccdef'
if val=='-1':
outColor='white'
if val=='ED0':
outColor='#FFFF00'
if val=='ED1':
outColor='#FCD37F'
if val=='ED2':
outColor='#FFAA00'
if val=='ED3':
outColor='#E60000'
if val=='ED4':
outColor='#730000'
ax.add_geometries([county], ccrs.PlateCarree(),facecolor=outColor, edgecolor='black',linewidth=0.10)
attribute_counter+=1
# Add Colormap
cmap = mpl.colors.ListedColormap(['#730000', '#E60000', '#FFAA00', '#FCD37F', '#FFFF00', 'white', '#8ccdef', '#00bfff', '#1d90ff', '#4169e1', '#0000ff'])
bounds=np.array([0,1,2,3,4,5,6,7,8,9,10,11],dtype='i')
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
cax = fig.add_axes([0.1, -0.035, 0.8, 0.03])
cbar=plt.colorbar(mpl.cm.ScalarMappable(cmap=cmap, norm=norm),cax=cax,boundaries=bounds,extend='neither',extendfrac='auto',ticks=bounds,spacing='uniform',orientation='horizontal')
cbar.ax.tick_params(labelsize=10)
cbar.set_label('Categories',size=15)
# Fix Tick Marks
loc = bounds + .5
cbar.set_ticks(loc)
labels=np.array(['ED4','ED3','ED2','ED1','ED0','None','EW0','EW1','EW2','EW3','EW4','EW5'],dtype='str')
cbar.set_ticklabels(labels)
# Add Titles
value_max="%6.1f" % np.nanmax(np.array(out_shapefile[inCode1]))
value_min="%6.1f" % np.nanmin(np.array(out_shapefile[inCode1]))
plt.suptitle("EDDI "+outElement+" Value For "+inDate,size=17,color='white',y=1.05)
#plt.annotate('MAX: '+str(value_max)+'\nMIN: '+str(value_min),xy=(0.045, -3.51), xycoords='axes fraction', fontsize=5,color='white',horizontalalignment='right', verticalalignment='bottom')
plt.annotate('Source: NOAA PSL\nMade By @jjrennie',xy=(1.045, -3.51), xycoords='axes fraction', fontsize=5,backgroundcolor='black',color='white',horizontalalignment='right', verticalalignment='bottom')
# Save Figure
plt.savefig(plot_directory+"/eddi-"+inLag+"_val_"+inDate+".png",bbox_inches='tight')
plt.clf()
plt.close()
#################################################
# PLOTTING CATEGORICAL VALUE
print("PLOTTING CATEGORICAL VALUE")
# Set Up Figure
fig= plt.figure(num=1, figsize=(8,5), dpi=dpi, facecolor='w', edgecolor='k')
ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.LambertConformal())
ax.set_facecolor('black')
ax.set_extent([minLon, maxLon, minLat, maxLat], crs=ccrs.Geodetic())
# Plot Data For Each County
attribute_counter=0
for county in shpreader.Reader(finalShapeFile).geometries():
val=out_shapefile.iloc[attribute_counter][inCode3]
if val != 'N/A':
if val=='-1':
outColor='#4a4a4a'
if val=='0':
outColor='#FFFF00'
if val=='1':
outColor='#FCD37F'
if val=='2':
outColor='#FFAA00'
if val=='3':
outColor='#E60000'
if val=='4':
outColor='#730000'
ax.add_geometries([county], ccrs.PlateCarree(),facecolor=outColor, edgecolor='black',linewidth=0.10)
attribute_counter+=1
# Add Colormap
cmap = mpl.colors.ListedColormap(['#FFFF00', '#FCD37F', '#FFAA00', '#E60000', '#730000'])
bounds=np.array([0,1,2,3,4,5],dtype='i')
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
cax = fig.add_axes([0.1, -0.035, 0.8, 0.03])
cbar=plt.colorbar(mpl.cm.ScalarMappable(cmap=cmap, norm=norm),cax=cax,boundaries=bounds,extend='neither',extendfrac='auto',ticks=bounds,spacing='uniform',orientation='horizontal')
cbar.ax.tick_params(labelsize=15)
cbar.set_label('USDM Drought Categories',size=15)
# Fix Tick Marks
loc = bounds + .5
cbar.set_ticks(loc)
labels=np.array(['D0','D1','D2','D3','D4','D5'],dtype='str')
cbar.set_ticklabels(labels)
# Add Titles
plt.suptitle("EDDI"+outElement+" USDM Category For "+inDate,size=17,color='white',y=1.05)
plt.annotate('Source: NOAA PSL\nMade By @jjrennie',xy=(1.045, -3.51), xycoords='axes fraction', fontsize=5,backgroundcolor='black',color='white',horizontalalignment='right', verticalalignment='bottom')
# Save Figure
plt.savefig(plot_directory+"/eddi-"+inLag+"_cat_"+inDate+".png",bbox_inches='tight')
plt.clf()
plt.close()
####################
# DONE
####################
print("DONE!")
end=time.time()
print(("Runtime: %8.1f seconds" % (end-start)))
sys.exit()
\ No newline at end of file
################################################################################
# Written By Jared Rennie
################################################################################
# Import Packages
import sys, time, datetime, os
import numpy as np
import numpy.ma as ma
import geopandas as gpd
import pandas as pd
import xarray as xr
import calendar
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
from scipy import stats
from PIL import Image
import warnings
warnings.filterwarnings("ignore")
# Define directories
main_directory="/store/sfcnet/jared/cpo_drought"
data_directory=main_directory+'/input_data/esi'
shapefile_directory=main_directory+'/results_shapefile'
csv_directory=main_directory+'/results_csv'
plot_directory=main_directory+'/results_plots'
#################################################
# BEGIN PROGRAM
start=time.time()
# Read in Arguments
if len(sys.argv) < 3:
sys.exit("USAGE: <Year-Month> <Month_Lag>\n example: ipython plot_esi.py 201608 01")
inDate=sys.argv[1]
inLag=sys.argv[2]
inCode1='esi_'+inLag+'-raw'
inCode2='esi_'+inLag+'-usdm'
if inLag=='01':
inWk='4WK'
outElement='(1-Month)'
if inLag=='03':
inWk='12WK'
outElement='(3-Month)'
# Get Day Info
lastDay=str(calendar.monthrange(int(inDate[0:4]), int(inDate[4:6]))[1])
fullDate=inDate+lastDay
dayStamps=np.array(['029','057','092','120','148','183','211','246','274','302','337','365'],dtype=str)
dayStamp=dayStamps[int(inDate[4:6])-1]
#################################################
# READ IN SHAPEFILE
input_shapefile=main_directory+'/input_shapefile/cb_2018_us_county_500k.shp'
print("READ IN SHAPEFILE: ",input_shapefile)
geo_shapefile = gpd.read_file(input_shapefile)
geo_shapefile['GEOID2_INT'] = np.array(geo_shapefile['GEOID'].values,dtype='i')
# Get Projection
projection=geo_shapefile.crs
#################################################
# READ IN GRID
num_lons=7200
num_lats=3000
min_lon=-180
min_lat=-60
cell_size=0.05
missing=-9999.
max_lon=180-cell_size
max_lat=90-cell_size
lons=np.linspace(min_lon,max_lon,num_lons)
lats=np.linspace(min_lat,max_lat,num_lats)
# Set Bounds (Gridding)
minLat = 20
maxLat = 55
minLon = -130
maxLon = -65
lats=np.flip(lats)
# Read In Data
input_grid=data_directory+'/DFPPM_'+str(inWk)+'_'+inDate[0:4]+dayStamp+'.tif'
print("READING IN DATA: "+str(input_grid))
im = Image.open(input_grid)
imarray = np.array(im)
# Ravel the Data to Match Lat/Lon Grid
gridlons, gridlats=np.meshgrid(lons,lats)
lats=gridlats.ravel()
lons=gridlons.ravel()
vals=imarray.ravel()
# Subset By Lat Bounds
latGrid=np.logical_and(lats >= minLat, lats <= maxLat)
lats=lats[latGrid]
lons=lons[latGrid]
vals=vals[latGrid]
# Subset By Lon Bounds
lonGrid=np.logical_and(lons >= minLon, lons <= maxLon)
lats=lats[lonGrid]
lons=lons[lonGrid]
vals=vals[lonGrid]
# Subset by NonMiss
nonMiss=np.where(vals!=missing)
lats=lats[nonMiss]
lons=lons[nonMiss]
vals=vals[nonMiss]
data_month = pd.DataFrame({'lat': lats, 'lon': lons, inCode1: vals})
# Convert table from Pandas to GeoPandas
geo_grid=gpd.GeoDataFrame(data_month,geometry=gpd.points_from_xy(data_month.lon, data_month.lat))
geo_grid['geometry']=geo_grid['geometry'].buffer(cell_size*4)
# Set Projection to same as Shapefile
geo_grid=geo_grid.set_crs(projection)
#################################################
# Spatial Join
print("SPATIAL JOIN")
# Perform the Join
result=gpd.sjoin(geo_shapefile,geo_grid,how="inner", op="intersects")
# Take the Aggregate Mean By County
result_mean=result.groupby(['GEOID'],as_index=False).mean()[['GEOID', inCode1]]
result_mean[inCode1]=result_mean[inCode1].round(decimals=2)
# Join Results with CONUS Counties, to reconcile missing counties during Spatial Join
conus_counties=pd.read_csv(main_directory+'/counties.csv')
conus_counties['GEOID']=conus_counties['GEOID'].astype(str).apply(lambda x: x.zfill(5))
merge_conus=conus_counties.merge(result_mean, on='GEOID', how='left')
result_mean=merge_conus
# Insert and Reformat Columns
result_mean[inCode1]=result_mean[inCode1].fillna(value=-9999.)
result_mean[inCode1]=result_mean[inCode1].round(decimals=1)
result_mean.insert(0, "Date", np.full(result_mean.shape[0],inDate), True)
result_mean.insert(3, inCode2, np.full(result_mean.shape[0],'-1'), True)
# Run through each county and Create USDM drought category.
for index, row in result_mean.iterrows():
val=row[inCode1]
if val==-9999.:
result_mean.at[index,inCode2]='N/A'
continue
category='-1'
if val<=-0.5:
category='0'
if val<=-0.8:
category='1'
if val<=-1.3:
category='2'
if val<=-1.6:
category='3'
if val<=-2.0:
category='4'
result_mean.at[index,inCode2]=category
# Create new shapfeile of Results
out_shapefile=geo_shapefile.merge(result_mean, on='GEOID')
# Save as new shapefile
finalShapeFile=shapefile_directory+"/esi-"+inLag+"_"+inDate+".shp"
print("OUTPUT TO: "+str(finalShapeFile))
out_shapefile.to_file(finalShapeFile)
# Fix SD FIPS so it can be 46102 in shapefile (to plot), but 46113 in csv (to match with health data)
fipsChange=result_mean.index[result_mean['GEOID'] == '46102'].tolist()[0]
result_mean.at[fipsChange,'GEOID']='46113'
# Save as new csv
finalTextFile=csv_directory+"/esi-"+inLag+"_"+inDate+".csv"
print("OUTPUT TO CSV: "+str(finalTextFile))
result_mean.to_csv(finalTextFile,index=False)
#################################################
# PLOTTING RAW VALUE
print("PLOTTING RAW VALUE")
unit='ESI Value'
color_map='BrBG'
bounds=np.linspace(-3.5,3.5,11)
dpi=300
# Set Bounds (Plotting)
minLat = 22
maxLat = 50
minLon = -120
maxLon = -73
# Set Up ColorBar
vmin=np.min(bounds)
vmax=np.max(bounds)
extend='both'
norm = mcolors.BoundaryNorm(boundaries=bounds, ncolors=256)
cm = plt.cm.ScalarMappable(cmap=color_map)
cm.set_array(np.array(out_shapefile[inCode1]))
cm.set_clim(vmin, vmax)
# Set Up Figure
fig= plt.figure(num=1, figsize=(8,5), dpi=dpi, facecolor='w', edgecolor='k')
ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.LambertConformal())
ax.set_extent([minLon, maxLon, minLat, maxLat], crs=ccrs.Geodetic())
plt.style.use('dark_background')
ax.set_facecolor('black')
# Plot Data For Each County
attribute_counter=0
for county in shpreader.Reader(finalShapeFile).geometries():
val=out_shapefile.iloc[attribute_counter][inCode1]
if val > -9999.:
ax.add_geometries([county], ccrs.PlateCarree(),facecolor=cm.to_rgba(val), edgecolor='black',linewidth=0.10)
attribute_counter+=1
# Add Colormap
cax = fig.add_axes([0.1, -0.035, 0.8, 0.03])
cbar=plt.colorbar(cm, cax=cax,boundaries=bounds,orientation='horizontal',extend=extend,spacing='uniform')
cbar.ax.tick_params(labelsize=15)
cbar.set_label(unit,size=15)
# Add Titles
value_max="%6.1f" % np.nanmax(np.array(out_shapefile[inCode1]))
value_min="%6.1f" % np.nanmin(np.array(out_shapefile[inCode1]))
plt.suptitle("ESI "+outElement+" Value For "+inDate,size=17,color='white',y=1.05)
#plt.annotate('MAX: '+str(value_max)+'\nMIN: '+str(value_min),xy=(0.045, -3.51), xycoords='axes fraction', fontsize=5,color='white',horizontalalignment='right', verticalalignment='bottom')
plt.annotate('Source: ClimateSERV\nMade By @jjrennie',xy=(1.045, -3.51), xycoords='axes fraction', fontsize=5,backgroundcolor='black',color='white',horizontalalignment='right', verticalalignment='bottom')
# Save Figure
plt.savefig(plot_directory+"/esi-"+inLag+"_val_"+inDate+".png",bbox_inches='tight')
plt.clf()
plt.close()
#################################################
# PLOTTING CATEGORICAL VALUE
print("PLOTTING CATEGORICAL VALUE")
# Set Up Figure
fig= plt.figure(num=1, figsize=(8,5), dpi=dpi, facecolor='w', edgecolor='k')
ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.LambertConformal())
ax.set_facecolor('black')
ax.set_extent([minLon, maxLon, minLat, maxLat], crs=ccrs.Geodetic())
# Plot Data For Each County
attribute_counter=0