path #= 0 length= 415.785621835 km
path #= 1 length= 527.959783762 km
path #= 2 length= 1432.08536051 km
path #= 3 length= 563.215908518 km
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# read data
data=pd.read_csv("chouisa.dat",sep='\s*',skiprows=1,header=None,index_col=0,parse_dates=True,names=("date","chouisa"),na_values="-")
#monthly value
data2=data.resample("M",how='mean')
#plot
plt.subplot(111)
data["chouisa"].plot(style='b-',label="Daily")
data2["chouisa"].plot(style='r-',label="Monthly",linewidth=2)
plt.legend(loc='best')
plt.savefig("choisa.png")
plt.show()
import numpy as np
from pandas import *
import matplotlib.pyplot as plt
ds = np.DataSource(None)
f=ds.open('http://www.data.kishou.go.jp/kaiyou/shindan/b_2/kuroshio_stream/kuro_slat.txt')
data=read_csv(f,sep=" ",names=("space","val"),parse_dates=True)
f.close()
data["val"].plot()
plt.xlabel("YEAR")
plt.ylabel("Kurosho latitude")
plt.savefig("pandas_example.png")
def get_gsm(date=[2011,1,1,0],dir="."):
"""
download jma gsm initial data from http://database.rish.kyoto-u.ac.jp/arch/jmadata/da
ta/gpv/original/
and
convert to grads file using gsm2bin by Dr. Shingo Shimizu
(http://shimizus.hustle.ne.jp/wiki/wiki.cgi?page=%A5%E2%A5%C7%A5%EB%B4%D8%CF%A2%A5%E1%A
5%E2)
usage:
get_gsm(date,dir)
input
date: [year,month,day,hour] hour must be multiple of 6
default [2011,1,1,0]
dir : directory to save data
"""
import os,urllib
#check
assert date[3] % 6 ==0, "Error, hour must be multiple of 6!!!"
#names from date
urldir,original,gradsbinname,gradsdate=name_maker(date)
# get original data
if not os.path.exists(dir+"/"+original):
if not os.path.exists(dir): os.mkdir(dir)
# Download the data
print 'Downloading data, please wait'
opener = urllib.urlopen(urldir+original)
open(dir+"/"+original, 'w').write(opener.read())
print 'Downloaded'
# convert to grads file
os.system("gsm2bin %(dir)s/%(original)s %(dir)s/%(gradsbinname)s" %locals())
# make ctlfile
make_ctl(dir,gradsbinname,gradsdate)
############################################
def name_maker(date):
"""
make filenames from date
"""
import datetime
time=datetime.datetime(*date)
cyear=time.strftime("%Y")
cmonth=time.strftime("%m")
cmonth2=(time.strftime("%b")).lower()
cday=time.strftime("%d")
chour=time.strftime("%H")
# URL of original data
urldir="http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/%(cyear)s/%(
cmonth)s/%(cday)s/" %locals()
# File name of original data
original="Z__C_RJTD_%(cyear)s%(cmonth)s%(cday)s%(chour)s0000_GSM_GPV_Rgl_FD0000_grib2.
bin" %locals()
# grads bin name after conversion
gradsbinname="gsm%(cyear)s%(cmonth)s%(cday)s%(chour)s.bin" %locals()
# date in grads format
gradsdate="%(chour)s:00Z%(cday)s%(cmonth2)s%(cyear)s" %locals()
return urldir,original,gradsbinname,gradsdate
###########################################
def make_ctl(dir,gradsbinname,gradsdate):
"""
make grads ctl file for GSM data
"""
#ctl file text
ctltxt="""
dset ^%(gradsbinname)s
options template
undef -999
xdef 720 LINEAR 0.0 0.5
ydef 361 LINEAR -90.0 0.5
zdef 17 LEVELS 1000 925 850 700 600 500 400 300 250 200 150 100 70 50
30 20 10
tdef 1 LINEAR %(gradsdate)s 6hr
vars 16
slp 0 0 sea level pressure
sp 0 0 surface pressure
su 0 0 surface westerly wind
sv 0 0 surface southerly wind
st 0 0 surface temperature
srh 0 0 surface relative humidity
lca 0 0 LowCloudAmount
mca 0 0 MidCloudAmount
hca 0 0 HighCloudAmount
tca 0 0 TotalCloudAmount
z 17 0 height
u 17 0 westerly wind
v 17 0 southerly wind
t 17 0 temperature
rh 17 0 relative humidity
w 17 0 p-velocity
endvars
"""
# write to file
gradsctlname=dir+"/"+gradsbinname.replace(".bin",".ctl")
f=open(gradsctlname,"w")
f.write(ctltxt %locals())
f.close()
#-----------------------------------------------------------------------
if __name__ == '__main__':
# get gsm data
get_gsm(date=[2011,1,1,6],dir="DATA")
# plot
from grads.galab import GaLab # for python interface for grads
ga = GaLab(Bin='grads',Window=False,Echo=True,Verb=0,Port=False)
ga.open("DATA/gsm2011010106.ctl")
script="""
set grads off
set gxout shaded
d slp/100
cbarn
draw title Sea level pressure (hpa)
printim gsm.png white
"""
ga(script)
import numpy as np
import os
def downloadfile(infile="remotefile",outfile="localfile"):
"""
download file from NEAR-GOOS data site with id and password
input
infile: remote file on the web
outfile: local file
set password and id below
"""
import urllib2
host="goos.kishou.go.jp:80"
realm="NEAR-GOOS"
usrname="xxxxx"
password="*****"
auth_handler = urllib2.HTTPBasicAuthHandler()
auth_handler.add_password(realm, host, usrname, password)
opener = urllib2.build_opener(auth_handler)
urllib2.install_opener(opener)
f=urllib2.urlopen(infile)
open(outfile, 'w').write(f.read())
#----------------------------------------------------------
def readMGDSST(dir=".",year=2010,month=1,day=1,lonrange=(0,360),latrange=(-90,90)):
"""
load MGDSST data from NEAR-GOOSE site
read data with fortran (f2py)
select region specified by lonrange and latrange
Usage:
lon,lat,sst= readMGDSST(dir=".",year=2010,month=1,day=1,lonrange=(0,360),latrange=(-90,90)
Input
dir: directory for data download
year,month,day: date of the data
lonrange,latrange = longitude and latitude range for data
output
lon, lat: longitude, latitude
sst: sst data
"""
# date
import datetime
t=datetime.datetime(year,month,day)
syear=t.strftime("%Y")
smon=t.strftime("%b")
sday=t.strftime("%d")
# filename
dataname="mgdsst."+smon+sday
filename=dir+"/"+dataname+"."+syear
filenamegz=filename+".gz"
print filename
urlname="http://goos.kishou.go.jp/rrtdb/usr/pub/JMA/mgdsst/"+syear+"/"+dataname+".gz"
# file exist?
if not os.path.exists(filename):
if not os.path.exists(filenamegz):
downloadfile(infile=urlname,outfile=filenamegz)
os.system("gunzip %(filenamegz)s" %locals())
# temp.bin is used for temporary file
if os.path.lexists("temp.bin"):
raise IOError("temporary file temp.bin somehow exist! Remove it!")
else:
os.system("ln -s %(filename)s temp.bin" %locals())
# data read through fortran (f2py)
from read_sst import read_sst
im=1440
jm=720
sst=np.zeros([jm,im])
sst,ierr=read_sst()
assert ierr==0, "data error in fortran"
lon=0.125+0.25*np.arange(im)
lat=-89.875+0.25*np.arange(jm)
os.system("rm temp.bin")
# range selection
xrange=np.logical_and(lon>=lonrange[0],lon<=lonrange[1])
yrange=np.logical_and(lat>=latrange[0],lat<=latrange[1])
sst=(sst[yrange,:])[:,xrange]
return lon[xrange],lat[yrange], sst
#-----------------------------------------------------------------------
if __name__ == '__main__':
import cdms2,MV2,vcs #cdat
# read data
lon,lat,sst=readMGDSST(dir="DATA",year=2011,month=1,day=1,lonrange=(0,360),latrange=(-90,90))
# hearafter plot with cdat
# First let's create a variable
var = MV2.masked_array(sst,mask=(sst<-990.),id="SST")
# Now create the axis
lons=cdms2.createAxis(lon)
lats=cdms2.createAxis(lat)
# And decorate it
lons.id = 'longitude'
lats.id = 'latitude'
lons.units = 'degrees_east'
lats.units = 'degrees_north'
# Just to make sure we will describe this axis as a longitude one
lons.designateLongitude()
lats.designateLatitude()
# Now let's apply these to the variable
var.setAxisList((lats,lons))
#var.info() to see info
#plot
v = vcs.init()
v.plot(var)
v.png("MGDSST.png")
del v
subroutine READ_SST(sst,ierr)
implicit none
integer :: year,mon,day
integer :: i,j,jj
integer,parameter :: im=1440
integer,parameter :: jm=720
!
real,intent(out):: sst(jm,im)
integer, intent(out):: ierr
!
integer :: isst(im,jm),ist
real :: undef = -999.
open(55,file="temp.bin",status='old',form='formatted',iostat=ierr,err=999)
read(55,*) year,mon,day
write(6,*) year,mon,day
do j=1,jm
read(55,10,iostat=ierr,err=999) (isst(i,j),i=1,im)
enddo
close(55)
10 format(1440i3)
do j=1,jm
do i=1,im
jj = jm -j +1
ist = isst(i,jj)
if(ist .eq. 888 .or. ist.eq.999) then
sst(j,i) = undef
else
sst(j,i) = real(ist)*0.1
endif
enddo
enddo
return
999 print *, "Read Error!"
end subroutine READ_SST
"reinit" *NOAA OISST V2 "sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/noaa.oisst.v2/sst.m nmean.nc" "set time jan2011" "mul 1 3 1 3 -yoffset -0.5" "set grads off" "color 0 30 2 -kind grainbow -gxout grfill" "d sst" "set strsiz 0.15" "draws 1) NOAA OISST V2" *NOAA Extended Reconstructed SST V3b "sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/noaa.ersst/sst.mnme an.nc" "set time jan2011" "mul 1 3 1 2 -yoffset -0.5 " "set grads off" "color 0 30 2 -kind grainbow -gxout grfill" "d sst.2" "set strsiz 0.15" "draws 2) NOAA Extended Reconstructed SST V3b" "cbarn 0.7 1 7.6 7.5" *diff "define sstre=re(sst.1,180,linear,0,2,89,linear,-88,2,ba)" "mul 1 3 1 1 -yoffset -0.5" "set grads off" "color -levs -2 -1.5 -1 -0.5 0.5 1.0 1.5 2 -kind blue->white->red -gxout grfill" "d sst.2-sstre" "set strsiz 0.15" "draws 2)-1)" "cbarn 0.7 1 1. 2." "printim sst_regrid_compare.png white"
def RadarToGrads(dir=".",gradsbinname="temp.bin",year=2010,mon=1,day=1,hour=0,mi
n=0,slon=118.,elon=150.,slat=20.,elat=48.):
"""
Import JMA rain radar data (every 10min) from http://database.rish.kyoto-u.ac
.jp/arch/jmadata/synthetic-original.html
Then convert data to grads file using tool by Dr. Shingo Shimizu (http://shim
izus.hustle.ne.jp/wiki/wiki.cgi?page=%A5%EC%A1%BC%A5%C0%B4%D8%CF%A2%A5%E1%A5%E2)
usage:
RadarToGrads(dir=".",gradsbinname="temp.bin",year=2010,mon=1,day=1,hour=0
,min=0,slon=118.,elon=150.,slat=20.,elat=48.):
input
dir: directory to import tar file
gradsbinname: file name for grads binary
year, mon, day, hour, min : time for data (every 10 min) in UTC
slon,elon Western,East bounday
slat,elat South,North bounday
output files
grads binary: gradsbinname
grads control: view_radar.ctl
"""
import os
import tarfile
import urllib
# filenames
syear="%04d" % year
smon="%02d" % mon
sday="%02d" % day
shour="%02d" % hour
smin="%02d" % min
date=syear+smon+sday+shour+smin
tarname="Z__C_RJTD_%(date)s00_RDR_JMAGPV__grib2.tar" %locals()
urldir="http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/jma-radar/synth
etic/original/%(syear)s/%(smon)s/%(sday)s/" %locals()
filename="Z__C_RJTD_%(date)s00_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_grib2.bin" %l
ocals()
# get tar file
if not os.path.exists(dir+"/"+tarname):
if not os.path.exists(dir): os.mkdir(dir)
# Download the data
print 'Downloading data, please wait'
opener = urllib.urlopen(urldir+tarname )
open(dir+"/"+tarname, 'w').write(opener.read())
print 'Downloaded'
# extract radar precipitation file
tarfile.open(dir+"/"+tarname, 'r').extract(filename)
# convert to grads file
os.system("jmaradar2bin %(filename)s %(gradsbinname)s %(slon)s %(elon)s %(sl
at)s %(elat)s" %locals())
os.system("rm %(filename)s" %locals())
####################################################################
def PlotByGrads(pngname):
"""
plot converted JMA radar rain data
"""
from grads.galab import GaLab # for python interface for grads
ga = GaLab(Bin='grads',Window=False,Echo=False,Verb=0,Port=False)
ga.open("view_radar.ctl")
script="""
set grads off
set mpdset hires
set gxout shaded
set clevs 0 1 5 10 20 30 50 80
d rr
cbarn
printim %(pngname)s white
"""
ga(script %locals())
#####################################################################
if __name__ == '__main__':
# radar file to grads
RadarToGrads(dir="DATA",gradsbinname="jmaradar.bin",year=2010,mon=10,day=20,
hour=8,min=00,slon=125.,elon=135.,slat=24.,elat=35.)
# plot by grads
PlotByGrads("jmaradar.png")
import numpy as np
# table parser from http://simbot.wordpress.com/2006/05/17/html-table-parser-using-python/
from table_parser import *
def read_FirstSnowFuji():
import urllib
"""
read first snow cover day of the FUJI mountain from KOFU local meteorological observatory
usage:
year,month,day=read_FirstSnowFuji()
input
none
output
year,month,day: first snow cover day (numpy array)
"""
f = urllib.urlopen('http://www.jma-net.go.jp/kofu/menu/siryo_hatsukansetsu.dat.html')
p = TableParser()
p.feed(f.read())
f.close()
data=p.doc # Get to the data
#get
year=[[],[],[]] # 3 column
month=[[],[],[]] # 3 column
day=[[],[],[]] # 3 column
for d in data[0][1:]:
for n,c in enumerate([0,4,8]):
if not d[c]=="":
year[n].append(int(d[c]))
mi,di=d[c+1].split("/")
month[n].append(int(mi))
day[n].append(int(di))
# to numpy array
ayear=np.array([],dtype="i4")
amonth=np.array([],dtype="i4")
aday=np.array([],dtype="i4")
for n in range(3):
ayear=np.append(ayear,year[n])
amonth=np.append(amonth,month[n])
aday=np.append(aday,day[n])
return ayear,amonth,aday
############################################
if __name__ == '__main__':
import matplotlib.pyplot as plt
import scikits.timeseries as ts
ayear,amonth,aday= read_FirstSnowFuji()
days=np.array([],dtype="i4")
for n in range(ayear.size):
td=ts.Date("D",year=ayear[n],month=amonth[n],day=aday[n])
days=np.append(days,td.day_of_year)
#plot
plt.plot(ayear,days,"b-")
mean=days.mean()
stdp=days.mean()+2*days.std()
stdm=days.mean()-2*days.std()
plt.plot([ayear[0],ayear[-1]],[mean,mean],"k-")
plt.plot([ayear[0],ayear[-1]],[stdp,stdp],"k--")
plt.plot([ayear[0],ayear[-1]],[stdm,stdm],"k--")
plt.xlim(ayear[0],ayear[-1])
plt.xlabel("year")
plt.ylabel("day of year")
plt.savefig("FujiFirstSnow.png")
plt.show()
import Nio
import numpy as np
def read_TRMM_3B42(date,lon1,lon2,lat1,lat2):
"""
read TRMM 3B42 3 hourly rainfall from NASA Mirador
usage:
lon,lat,sst=read__TRMM_3B42(date,lon1,lon2,lat1,lat2)
input
date (year,month,day,hour)
lon1,lon2 Western,East bounday
lat1,lat2 South,West bounday
output
lon longitude
lat latitude
rainfall data (mm/hour)
"""
#check
assert date[3] % 3 ==0, "Error, hour must be multiple of 3!!!"
# open OPENDAP
fname=name_maker(date[0],date[1],date[2],date[3])
f = Nio.open_file(fname)
## variable name check
#variables = f.variables.keys()
#print f
#print "variables",variables
# read data
temp=f.variables["DATA_GRANULE.PlanetaryGrid.precipitation"]
#grid
# grid
lon=-179.875+np.arange(temp.shape[1])*0.25
lat=-49.875 +np.arange(temp.shape[2])*0.25
# selection
xrange=np.where(np.logical_and(lon>=lon1,lon<=lon2))
yrange=np.where(np.logical_and(lat>=lat1,lat<=lat2))
#precipitation
x=xrange[0]
y=yrange[0]
precip=(temp[0,x[0]:x[-1]+1,y[0]:y[-1]+1].T).reshape((len(y),len(x)))
# closing
f.close()
return lon[xrange],lat[yrange],precip
############################################
def name_maker(year,month,day,hour):
import scikits.timeseries as ts
#
t=ts.Date('H', year=year,month=month,day=day,hour=hour)
if t.hour==0:
tc=t-1
else:
tc=t
cyear = "%04d" % tc.year
cyear2= "%02d" % (t.year-2000)
cmonth2= "%02d" % t.month
cday2= "%02d" % t.day
shour = str(t.hour)
yd3= "%03d" % tc.day_of_year
#
fname="http://disc2.nascom.nasa.gov/opendap/TRMM_L3//TRMM_3B42//%(cyear)s/%(yd3)s/3B42.%(cyear2)s%(cmonth2)s%(cday2)s.%(shour)s.6A.HDF.Z" %locals()
return fname
############################################
if __name__ == '__main__':
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
# example plot
# range
lon1=120.
lon2=145.
lat1=20.
lat2=40.
#date
year=2009
month=7
day=27
hour=0
lon,lat,precip=read_TRMM_3B42(date=(year,month,day,hour),lon1=lon1,lon2=lon2,lat1=lat1,lat2=lat2)
#plot
vmax=int(precip.max())
vmin=int(precip.min())
fig=plt.figure()
ax = fig.add_axes([0.1,0.1,0.8,0.8])
map = Basemap(projection='merc',llcrnrlat=lat1,urcrnrlat=lat2,\
llcrnrlon=lon1,urcrnrlon=lon2,resolution='i')
lons,lats=np.meshgrid(lon[:],lat[:])
x,y=map(lons,lats)
map.drawcoastlines()
map.drawmapboundary()
map.pcolor(x,y,precip,vmin=vmin,vmax=vmax,cmap=plt.cm.Blues)
ti=str(year)+"/"+str(month)+"/"+str(day)+":"+str(hour)+"Z"
plt.title(ti)
# colorbar
pos = ax.get_position()
l, b, w, h = pos.bounds
cax = plt.axes([l+w+0.02, b, 0.03, h]) # setup colorbar axes.
plt.colorbar(cax=cax) # draw colorbar
plt.savefig("TRMM_3B42.png")
plt.show()
import numpy as np
def read_MODIS_EORC(file,lon1,lon2,lat1,lat2):
"""
read MODIS_SST from EORC
usage:
lon,lat,sst=read_MODIS_EORC(file,lon1,lon2,lat1,lat2)
input
file filename
lon1,lon2 Western,East bounday
lat1,lat2 South,West bounday
output
lon longitude
lat latitude
SST sst data (masked array) in C
"""
# open file
f=open(file, 'r')
# read header
pixel=f.read(6); print "pixel=",pixel
line=f.read(6); print "line=",line
nwlon=f.read(8); print "NW lon.=",nwlon
nwlat=f.read(8); print "NW lat.=",nwlat
res=f.read(8); print "resolution=",res
slope=f.read(9); print "slope=",slope
offset=f.read(9); print "ofset=",offset
dummy=f.read(1)
parameter=f.read(8); print "parameter name=",parameter
dummy=f.read(1)
filename=f.read(45); print "file name=",filename
f.close()
# read data
f=open(file,'r')
f.seek(int(pixel)*2) #skip header
data=np.fromfile(f,dtype=">u2")
data=data.reshape((int(line),int(pixel)))
msk=np.logical_or(data==0,data>=65534) # 0 cloud, 65534 out of obs., 65535 land
data=np.ma.masked_array(data[::-1],mask=msk[::-1,:])
data=data*float(slope)+float(offset)-273.15 # to degree C
f.close()
# grid
lon=float(nwlon)+np.arange(int(pixel))*float(res)
lat=float(nwlat)-np.arange(int(line))*float(res)
lat=lat[::-1]
# selection
xrange=np.logical_and(lon>=lon1,lon<=lon2)
yrange=np.logical_and(lat>=lat1,lat<=lat2)
sst=(data[yrange,:])[:,xrange]
return lon[xrange],lat[yrange], sst
############################################
if __name__ == '__main__':
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
# example plot
#read data
filename="A2GL11008020431OD1_OSTAQ_01000_01000_sst_06"
lon1=135.
lon2=138.
lat1=32.
lat2=35.
lon,lat,sst=read_MODIS_EORC(filename,lon1,lon2,lat1,lat2)
#plot
tmax=int(sst.max())
tmin=int(sst.min())
fig=plt.figure()
ax = fig.add_axes([0.1,0.1,0.8,0.8])
map = Basemap(projection='merc',llcrnrlat=lat1,urcrnrlat=lat2,\
llcrnrlon=lon1,urcrnrlon=lon2,resolution='i')
lons,lats=np.meshgrid(lon[:],lat[:])
x,y=map(lons,lats)
map.drawcoastlines()
map.drawmapboundary()
map.pcolor(x,y,sst,vmin=tmin,vmax=tmax)
# colorbar
pos = ax.get_position()
l, b, w, h = pos.bounds
cax = plt.axes([l+w+0.02, b, 0.03, h]) # setup colorbar axes.
plt.colorbar(cax=cax) # draw colorbar
plt.savefig("MODIS_EORC_SST.png")
plt.show()
import numpy as np
import matplotlib.pyplot as plt
def read_CEReS_SST(dir,file,lon1,lon2,lat1,lat2):
"""
read CERes SST
usage:
lon,lat,sst=read_CEReS_SST(dir,file,lon1,lon2,lat1,lat2)
input
dir working dir
file filename
gzipped file is fine.
URL also works
lon1,lon2 Western,East bounday
lat1,lat2 South,West bounday
output
lon longitude
lat latitude
SST sst data (masked array) in C
"""
# read data
ds=np.DataSource(dir)
f=ds.open(file, 'r')
f.seek(80) # skip header
ssttemp=np.frombuffer(f.read(),dtype=">i2",count=6378*5562)
f.close()
# as masked array
sst=ssttemp.reshape((5562,6378))*0.1
sst=np.ma.masked_array(sst[::-1,:],mask=(sst[::-1,:]==0 ))
# grid
lon=100.+np.arange(6378)*0.01097869
lat=9.97971+np.arange(5562)*0.00899322
# selection
xrange=np.logical_and(lon>=lon1,lon<=lon2)
yrange=np.logical_and(lat>=lat1,lat<=lat2)
sst=(sst[yrange,:])[:,xrange]-273.15 # Kelvin to C
return lon[xrange],lat[yrange], sst
############################################
if __name__ == '__main__':
from mpl_toolkits.basemap import Basemap
# example plot
#read data
filename="n1810080316.sst.gi"
lon1=130.
lon2=140.
lat1=30.
lat2=36.
lon,lat,sst=read_CEReS_SST(".",filename,lon1,lon2,lat1,lat2)
#plot
tmax=30.
tmin=20.
fig=plt.figure()
ax = fig.add_axes([0.1,0.1,0.8,0.8])
map = Basemap(projection='merc',llcrnrlat=lat1,urcrnrlat=lat2,\
llcrnrlon=lon1,urcrnrlon=lon2,resolution='i')
lons,lats=np.meshgrid(lon[:],lat[:])
x,y=map(lons,lats)
map.drawcoastlines()
map.drawmapboundary()
map.pcolor(x,y,sst,vmin=tmin,vmax=tmax)
# colorbar
pos = ax.get_position()
l, b, w, h = pos.bounds
cax = plt.axes([l+w+0.02, b, 0.03, h]) # setup colorbar axes.
plt.colorbar(cax=cax) # draw colorbar
plt.savefig("CEReS_SST.png")
plt.show()
