Commit 30615d57 authored by SemjonSchimanke's avatar SemjonSchimanke
Browse files

Initial check-in of maps.

parent 08dccce2
This repository is a collection of tools for the PRECISE service running at SMHI.
Structure:
* maps: plotting routines for various parameters as verification
for the operational production
- parameters: T2m, ...
* ERA-diff: plotting routines showing the difference between UERRA/PRECISE and ERA-interim/ERA5
"""
A collection of useful colormaps.
"""
def red_white_blue_cmap(number):
"""
Create a color map from red to blue
with white in the middle. You have to define the number of
color steps you want.
Returns a cmap-class.
Example for usage of the returned cmap-class:
m.contourf(x,y,data,cmap=my_cmap)
"""
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
if not isinstance(number, int):
print "Please give an integer between 10 and 100."
return None
if number < 10 or number > 100:
print "Please give an integer between 10 and 100."
return None
jetcmap = plt.cm.get_cmap("jet", number) #generate a jet map with "number" values
jet_vals = jetcmap(np.arange(number)) #extract those values as an array
# Make the colors around middle white
nn=number/2
if number < 21:
jet_vals[nn-1] = [1, 1, 1, 1] #change a value
jet_vals[nn] = [1, 1, 1, 1]
else:
jet_vals[nn-2] = [1, 1, 1, 1] #change a value
jet_vals[nn-1] = [1, 1, 1, 1]
jet_vals[nn] = [1, 1, 1, 1]
jet_vals[nn+1] = [1, 1, 1, 1]
return mpl.colors.LinearSegmentedColormap.from_list("newjet", jet_vals)
def white_blue_red_cmap(number):
"""
Create a color map from blue to red
with white at the lower end. You have to define the number of
color steps you want.
Returns a cmap-class.
"""
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
if not isinstance(number, int):
print "Please give an integer between 10 and 100."
return None
if number < 10 or number > 100:
print "Please give an integer between 10 and 100."
return None
jetcmap = plt.cm.get_cmap("jet", number) #generate a jet map with "number" values
jet_vals = jetcmap(np.arange(number)) #extract those values as an array
# Change first colors to be white
jet_vals[0] = [1, 1, 1, 1] #change a value
jet_vals[1] = [1, 1, 1, 1]
return mpl.colors.LinearSegmentedColormap.from_list("newjet", jet_vals)
# blue to red
# This can be used with "colors=RGB_tuples"
blue_to_red_tuples_18=[(1.0, 1.0, 1.0),
(0.1, 0.1, 1.0),
(0.2, 0.2, 1.0),
(0.3, 0.3, 1.0),
(0.4, 0.4, 1.0),
(0.5, 0.5, 1.0),
(0.6, 0.6, 1.0),
(0.7, 0.7, 1.0),
(0.8, 0.8, 1.0),
(1.0, 0.8, 0.8),
(1.0, 0.7, 0.7),
(1.0, 0.6, 0.6),
(1.0, 0.5, 0.5),
(1.0, 0.4, 0.4),
(1.0, 0.3, 0.3),
(1.0, 0.2, 0.2),
(1.0, 0.1, 0.1),
(0.5, 1.0, 0.0)]
blue_to_red_tuples_24=[(0.0, 0.0 , 1.0),
(0.1, 0.1, 1.0),
(0.0, 0.1667, 1.0),
(0.2, 0.2, 1.0),
(0.3, 0.5 , 1.0),
(0.3, 0.7 , 1.0),
(0.3, 0.8333, 1.0),
(0.4, 0.95 , 1.0),
(0.3, 1.0 , 0.9300),
(0.0, 1.0 , 0.7125),
(0.0, 1.0 , 0.5125),
(0.5, 1.0 , 0.5),
(0.9, 1.0 , 0.2),
(1.0, 1.0 , 0.0),
(1.0, 0.8330, 0.0),
(1.0, 0.6667, 0.0),
(1.0, 0.3167, 0.0),
(1.0, 0.2, 0.2),
(1.0, 0.0 , 0.0),
(0.8, 0.1 , 0.5),
(0.541, 0.169, 0.886),
(0.282, 0.239, 0.545),
(0.4, 0.0 , 0.0),
(0.,0.,0.)]
white_blue_red_tuples=[(1.0, 1.0 , 1.0),
(0.1, 0.1, 1.0),
(0.0, 0.1667, 1.0),
(0.2, 0.2, 1.0),
(0.3, 0.5 , 1.0),
(0.3, 0.7 , 1.0),
(0.3, 0.8333, 1.0),
(0.4, 0.95 , 1.0),
(0.3, 1.0 , 0.9300),
(0.0, 1.0 , 0.7125),
(0.0, 1.0 , 0.5125),
(0.5, 1.0 , 0.5),
(0.9, 1.0 , 0.2),
(1.0, 1.0 , 0.0),
(1.0, 0.8330, 0.0),
(1.0, 0.6667, 0.0),
(1.0, 0.3167, 0.0),
(1.0, 0.2, 0.2),
(1.0, 0.0 , 0.0),
(0.8, 0.1 , 0.5),
(0.541, 0.169, 0.886),
(0.282, 0.239, 0.545),
(0.4, 0.0 , 0.0),
(0.,0.,0.)]
#!/usr/bin/python
#SBATCH -t 0-1
# from Semjon S.:
# program to loop over files and forecast lengths fr further
# handing off to plotting routines
import datetime
import subprocess
#firstdate=datetime.date(2014,12,2)
#lastdate =datetime.date(2015,12,31)
datelist=[2016020112]
fc_hours=['00', '01', '02', '03', '04', '05', '06', '09', '12', '18']
for date in datelist:
for hour in fc_hours:
subprocess.check_call(["python", "plot2png_diff_uerra.py", str(date),
hour])
#python plot2png_diff_uerra.py 2016020112 03
#!/bin/bash
#SBATCH --qos=normal
#SBATCH --export=NONE
#SBATCH --error=job_maps.log
#SBATCH --job-name=UerraMps
#SBATCH --output=job_maps.log
#SBATCH --workdir=/home/ms/se/smu/maps
#SBATCH --time=00:30:00
module load python/2.7.12-01
dat=$1
cd $SCRATCH/maps/
# plot maps for the long-run cycles only, i.e. 00UTC and 12UTC
for hh in '00' '12';
do
dats=$1$hh
echo 'plotting maps for ' $dats
for LL in '00' '03' '06' '09' '12' '15' '18' '21' '24' '27' '30' ;
do
python $HOME/maps/plot2png_uerra_surflevs.py $dats $LL
python $HOME/maps/plot2png_uerra_upperlevs.py $dats $LL
done
done
# now make sure the maps are placed in the right sub-directories
for vars in '10mWind' '300hPaWind' 'CClCov' 'gp500_mslp' 'HClCov' 'LClCov' 'MClCov' \
'MSLP-GP500-PPN' 'MSLP-PPN' 'precip' 'Rh2m' 'rh850' 'T2m' 't850' 'TClCov' ;
do
mkdir -p $dat/$vars
mv *$dat*$vars.png $dat/$vars/.
done
## Finally might want to tar and gzip the directory and maps for the day...
## tar -czf 20171104.tar.gz 20171104
#
tar -czf $HOME/maps/$dat.tar.gz $dat
cd $HOME/maps
#!/bin/bash
#SBATCH --qos=normal
#SBATCH --export=NONE
#SBATCH --error=job_maps.log
#SBATCH --job-name=UerraMps
#SBATCH --output=job_maps.log
#SBATCH --workdir=/home/ms/se/smu/maps
#SBATCH --time=00:30:00
module load python/2.7.12-01
dat=$1
cd $SCRATCH/maps/
# plot maps for the long-run cycles only, i.e. 00UTC and 12UTC
for hh in '00' '12';
do
dats=$1$hh
echo 'plotting maps for ' $dats
for LL in '00' '03' '06' '09' '12' '15' '18' '21' '24' '27' '30' ;
do
python $HOME/maps/plot2png_uerra_surfflx_fp.py $dats $LL
done
done
# now make sure the maps are placed in the right sub-directories
for vars in '000-111' '000-112' '000-132' '000-122'
do
mkdir -p $dat/$vars
mv *$dat*$vars.png $dat/$vars/.
done
## Finally might want to tar and gzip the directory and maps for the day...
## tar -czf 20171104.tar.gz 20171104
#
tar -czf $HOME/maps/$dat.fluxes.tar.gz $dat
cd $HOME/maps
#!/usr/bin/env python
# Import modules
import os
import glob
import sys
import pygrib
from pylab import *
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
import color_map
# used to construct filename to requested date, these need to be strings
# cdate = 'yyyymmddhh' & 'll' = length of forecast (2 digits)
# e.g. call with
# python plot2png_diff_uerra.py 2016010212 01
print sys.argv[1],sys.argv[2]
cdate = sys.argv[1]
ll = sys.argv[2]
def getgrbfield(name1,typeoflev1,level1):
selgrbs = grbindx(name=name1,typeOfLevel=typeoflev1,level=level1)
print selgrbs
datafield=selgrbs[0].values
return datafield, selgrbs
def plotanydiff(x,y,quant,titletext,levels,figno):
"""
Plotting function for the parameters
T2m, MSLP, RH2m, total precipitation.
Total Cloud, high-, mid- low clouds, convective clouds
The DIFFERENCES between 2 experiments is plotted.
You need to adjust the filenames to suit your needs.
"""
plt.figure(figno)
plt.clf()
rwb_cmap=color_map.red_white_blue_cmap(len(levels)-1)
m.contourf(x,y,quant, levels, cmap=rwb_cmap, extend='both')
m.drawcoastlines(color='gray')
m.drawcountries(color='lightgray')
m.drawmapboundary(color='lightgray')
xlabel('Lon')
ylabel('Lat')
title(titletext+cdate+' '+ll+'h')
plt.colorbar()
plt.savefig('./fc_fp_diff'+cdate+'+0'+ll+titletext+'.png')
# plt.show()
#
# # # ===============================================================
# Settings of variable name and level. This output will be also generated below.
name1 = ('Total precipitation')
level1 = 0
typeoflev1 = ('heightAboveGround')
name2 = ('u-component of wind','v-component of wind')
#name2 = ('U component of wind','V component of wind') # That's the name on Bi.
level2 = 10
typeoflev2 = ('heightAboveGround')
name3 = ('Pressure')
level3 = 0
typeoflev3 = ('heightAboveSea')
name4 = ('Temperature')
level4 = 2
typeoflev4 = ('heightAboveGround')
name5 = ('Relative humidity')
level5 = 2
typeoflev5 = ('heightAboveGround')
name6 = ('Total cloud cover')
level6 = 0
typeoflev6 = ('heightAboveGround')
name7 = ('High cloud cover')
level7 = 0
typeoflev7 = ('heightAboveGround')
name8 = ('Medium cloud cover')
level8 = 0
typeoflev8 = ('heightAboveGround')
name9 = ('Low cloud cover')
level9 = 0
typeoflev9 = ('heightAboveGround')
name10 = ('Convective cloud cover')
level10 = 0
typeoflev10 = ('heightAboveGround')
# Name of grib file, just for 10m data
inFile1 = '/scratch/ms/se/smu/tmp/uerra_ald_2016/fc'+cdate+'+0'+ll+'grib_fp'
inFile2 = '/scratch/ms/se/smu/tmp/uerra_ald_tke/fc'+cdate+'+0'+ll+'grib_fp'
# the filepath on Bi
#inFile1 = '/nobackup/fouo5/sm_saefa/UERRA/data/uerra_ald_2016/fc'+cdate+'+0'+ll+'grib_fp'
#inFile2 = '/nobackup/fouo5/sm_saefa/UERRA/data/uerra_ald_tke/fc'+cdate+'+0'+ll+'grib_fp'
grbs1 = pygrib.open(inFile1)
# Create grib index instance for faster searching
grbindx = pygrib.index(inFile1,'name','typeOfLevel','level')
# Get the data
precip_1,selgrbs1 = getgrbfield(name1,typeoflev1,level1)
u10_1,selgrbs1 = getgrbfield(name2[0],typeoflev2,level2)
v10_1,selgrbs1 = getgrbfield(name2[1],typeoflev2,level2)
mslp_1,selgrbs1 = getgrbfield(name3,typeoflev3,level3)
t2m_1,selgrbs1 = getgrbfield(name4,typeoflev4,level4)
rh2m_1,selgrbs1 = getgrbfield(name5,typeoflev5,level5)
tclcov_1,selgrbs1 = getgrbfield(name6,typeoflev6,level6)
hclcov_1,selgrbs1 = getgrbfield(name7,typeoflev7,level7)
mclcov_1,selgrbs1 = getgrbfield(name8,typeoflev8,level8)
lclcov_1,selgrbs1 = getgrbfield(name9,typeoflev9,level9)
cclcov_1,selgrbs1 = getgrbfield(name10,typeoflev10,level10)
# close grib index
grbindx.close()
windspeed_1=np.sqrt(u10_1**2 + v10_1**2)
grbs2 = pygrib.open(inFile2)
# Create grib index instance for faster searching
grbindx = pygrib.index(inFile2,'name','typeOfLevel','level')
# Get the data
precip_2,selgrbs2 = getgrbfield(name1,typeoflev1,level1)
u10_2,selgrbs2 = getgrbfield(name2[0],typeoflev2,level2)
v10_2,selgrbs2 = getgrbfield(name2[1],typeoflev2,level2)
mslp_2,selgrbs2 = getgrbfield(name3,typeoflev3,level3)
t2m_2,selgrbs2 = getgrbfield(name4,typeoflev4,level4)
rh2m_2,selgrbs2 = getgrbfield(name5,typeoflev5,level5)
tclcov_2,selgrbs2 = getgrbfield(name6,typeoflev6,level6)
hclcov_2,selgrbs2 = getgrbfield(name7,typeoflev7,level7)
mclcov_2,selgrbs2 = getgrbfield(name8,typeoflev8,level8)
lclcov_2,selgrbs2 = getgrbfield(name9,typeoflev9,level9)
cclcov_2,selgrbs2 = getgrbfield(name10,typeoflev10,level10)
# close grib index
grbindx.close()
windspeed_2=np.sqrt(u10_2**2 + v10_2**2)
# create a shortcut for the selected data
grb1 = selgrbs1[0]
grb2 = selgrbs2[0]
# Get the data and the lat/lon values
data1=grb1.values
data2=grb2.values
#differences for the various variables
precipdiff=precip_1-precip_2
u10diff=u10_1-u10_2
v10diff=v10_1-v10_2
windspeeddiff=windspeed_1-windspeed_2
mslpdiff=mslp_1-mslp_2
t2mdiff=t2m_1-t2m_2
rh2mdiff=rh2m_1-rh2m_2
# and differences in clouds (total,high,mid,low,convective)
tclcvdiff=tclcov_1-tclcov_2
hclcvdiff=hclcov_1-hclcov_2
mclcvdiff=mclcov_1-mclcov_2
lclcvdiff=lclcov_1-lclcov_2
cclcvdiff=cclcov_1-cclcov_2
#data=grb1.values-grb2.values
data=data1-data2
print 'shape/min/max data %s %6.2f %6.2f'%(str(data1.shape),data1.min(),data1.max())
lats, lons = grb1.latlons() # returns lat/lon values on grid.
print str('min/max of %d lats on %s grid %4.2f %4.2f' % (grb1['Nj'], grb1['typeOfGrid'],lats.min(),lats.max()))
print str('min/max of %d lons on %s grid %4.2f %4.2f' % (grb1['Ni'], grb1['typeOfGrid'],lons.min(),lons.max()))
# Create map, lowerleftcorner longitude and latitude, upperrightcorner ...
llcrnrlon = 0
llcrnrlat = 54
urcrnrlon = 30
urcrnrlat = 78
#
lon_0 = 10
lat_0 = 49
width = 8000000
height= 8000000
rsphere = (grb1.projparams['a'], grb1.projparams['b'])
projection = grb1.projparams['proj']
# project the data on the above map
m = Basemap(lon_0=lon_0,lat_0=lat_0,width=width,height=height,rsphere=rsphere,
resolution='l',projection='lcc')
x,y = m(lons, lats)
## Plot of field
## Put interactive mode on for plotting
#plt.ion()
figno=0
# Start plot festival
levels=np.linspace(-5,5,21)
plotanydiff(x,y,t2mdiff,'T2m',levels,figno+1)
levels=np.linspace(-0.2,0.2,21)
plotanydiff(x,y,rh2mdiff,'Rh2m',levels,figno+2)
levels=np.linspace(-10,10,21)
plotanydiff(x,y,windspeeddiff,'Windspeed',levels,figno+3)
levels=np.linspace(-0.5,0.5,21)
plotanydiff(x,y,precipdiff,'Precipitation',levels,figno+4)
levels=np.linspace(-20,20,21)
plotanydiff(x,y,mslpdiff,'MSLP',levels,figno+5)
levels=np.linspace(-0.5,0.5,21)
plotanydiff(x,y,tclcvdiff,'TClCov',levels,figno+6)
plotanydiff(x,y,hclcvdiff,'HClCov',levels,figno+7)
plotanydiff(x,y,mclcvdiff,'MClCov',levels,figno+8)
plotanydiff(x,y,lclcvdiff,'LClCov',levels,figno+9)
plotanydiff(x,y,cclcvdiff,'CClCov',levels,figno+10)
###levels=np.linspace(-10,10,21)
###plotwinddiff(x,y,u10diff,v10diff,levels,figno+3)
#plotprecipdiff(x,y,precipdiff,figno+4)
#plotmslpdiff(x,y,mslpdiff,figno+5)
#raw_input("Press enter to exit")
#!/usr/bin/env python
# Import modules
import os
import glob
import sys
import pygrib
from pylab import *
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm
import numpy as np
#import pylab
##yr,mo,da,hr= (2016,1,01,12)
# used to construct filename to requested date, these need to be strings
# cdate = 'yyyymmddhh' & 'll' = length of forecast (2 digits)
print sys.argv[1],sys.argv[2]
cdate = sys.argv[1]
ll = sys.argv[2]
def getgrbfield(name1,typeoflev1,level1):
selgrbs = grbindx(name=name1,typeOfLevel=typeoflev1,level=level1)
print selgrbs
datafield=selgrbs[0].values
return datafield, selgrbs
def plotany(x,y,quant,titletext,figno):
plt.figure(figno)
plt.clf()
if titletext == 'T2m':
clevs=np.arange(-50,50,2)
cs=m.contourf(x,y,quant,clevs,cmap=plt.cm.RdBu_r)
plt.colorbar(cs)
elif titletext == 'Rh2m':
clevs=np.arange(0.05,1.1,0.02)
cs=m.contourf(x,y,quant,clevs,cmap=plt.cm.gist_earth_r)
plt.colorbar(cs)
elif titletext == 'TClCov':
clevs=np.arange(0.05,1.1,0.02)
cs=m.contourf(x,y,quant,clevs,cmap=plt.cm.bone_r)
plt.colorbar(cs)
elif titletext == 'HClCov':
clevs=np.arange(0.05,1.1,0.02)
cs=m.contourf(x,y,quant,clevs,cmap=plt.cm.Blues)
plt.colorbar(cs)
elif titletext == 'MClCov':
clevs=np.arange(0.05,1.1,0.02)
cs=m.contourf(x,y,quant,clevs,cmap=plt.cm.Greens)
plt.colorbar(cs)
elif titletext == 'LClCov':
clevs=np.arange(0.05,1.1,0.02)
cs=m.contourf(x,y,quant,clevs,cmap=plt.cm.YlOrRd)
plt.colorbar(cs)
elif titletext == 'CClCov':
clevs=np.arange(0.05,1.1,0.02)
cs=m.contourf(x,y,quant,clevs,cmap=plt.cm.hot_r)
plt.colorbar(cs)
else:
m.contourf(x,y,quant,30)
plt.colorbar()
# m.contourf(x,y,quant,30)
m.drawcoastlines(color='gray')
m.drawcountries(color='lightgray')
m.drawmapboundary(color='lightgray')
xlabel('Lon')
ylabel('Lat')
title(titletext+cdate+' '+ll+'h')
pngfile='fc_fp_'+cdate+'+0'+ll+'_'+titletext+'.png'
# plt.colorbar()
plt.savefig(pngfile)
plt.show()
def plotprecip(x,y,precip,figno):
plt.figure(figno)
precip[precip<0.01] = NaN
#clevs = [0,0.01,1.0,2.5,5.0,7.5,10,15,20,30,40,50,70,80,90,100.,150.]
#clevs = [0,0.05,0.1,0.2,0.5,0.75,1.0,1.5,2.,3.,4.,5.,7.5,10,15,20,25,30,40,50,70,80,90,100.,150.]
#cs=m.contourf(x,y,precip,clevs,cmap=cm.s3pcpn)
clevs = [0.1,0.2,0.5,0.75,1.0,1.5,2.,3.,4.,5.,7.5,10,15,20,25,30,40,50,70,80,90,100.,150.]
cs=m.contourf(x,y,precip,clevs,cmap=plt.cm.Spectral_r)
m.drawcoastlines(color='gray')
m.drawcountries(color='lightgray')
m.drawmapboundary(color='lightgray')
xlabel('Lon')
ylabel('Lat')
title('Total precipitation'+cdate+' '+ll+'h')
plt.colorbar(cs)
plt.savefig('./fc_fp_'+cdate+'+0'+ll+'_precip.png')
plt.show()
def plotwind(x,y,u,v,figno):
ws= np.sqrt(u**2 + v**2)
skiparrow=10
plt.figure(figno)
clevs = np.arange(5,60)
cs=m.contourf(x,y,ws,clevs</