Hi All,
I have monthly netcdf files containing hourly temperature data.
I would like to loop through all of the hours and calculate the max for
each hour (00 - 23) and then make a plot.
I have found a way for this to work for a couple of the hours (see below)
but it requires a lot of set up to do each hour and I am sure there must be
a simpler way using loops? I also need to find a simpler way as I would
like to eventually do daily for 1st jan to 31st Dec
Any feedback will be greatly appreciated
from netCDF4 import Dataset
import numpy as N
import matplotlib.pyplot as plt
from numpy import ma as ma
from mpl_toolkits.basemap import Basemap
from netcdftime import utime
from datetime import datetime
import os
from matplotlib.collections import LineCollection
shapefile1="---"
OutputFolder=r"/---/"
fileforlatlon=Dataset("---.nc", 'r+', 'NETCDF4')
LAT=fileforlatlon.variables['latitude'][:]
LON=fileforlatlon.variables['longitude'][:]
#Set up basemap using mercator projection
http://matplotlib.sourceforge.net/basemap/doc/html/users/merc.html
map =
Basemap(projection='merc',llcrnrlat=-40,urcrnrlat=-33,llcrnrlon=139.0,urcrnrlon=151.0,lat_ts=0,resolution='i')
x,y=map(*N.meshgrid(LON,LAT))
map.readshapefile(shapefile1, '-REGIONS')
ncvariablename='T_SFC'
MainFolder=r"/---/"
ticks=[-5,0,5,10,15,20,25,30,35,40,45,50]
Title='Surface Temperature (degrees celsius)'
cmap=plt.cm.jet
all_variabledata=[]
time00=[]
time12=[]
for (path, dirs, files) in os.walk(MainFolder):
for dir in dirs:
print "the dirs are:", dir
path=path+'/'
for ncfile in files:
fileext=ncvariablename+'.nc'
if ncfile.endswith(fileext):
print "dealing with ncfiles:", path+ncfile
ncfile=os.path.join(path,ncfile)
ncfile=Dataset(ncfile, 'r+', 'NETCDF4')
TIME=ncfile.variables['time'][1::]
#variable=ncfile.variables[ncvariablename][:,:,:]
TSFC00=ncfile.variables[ncvariablename][00:01:,:,:]
TSFC12=ncfile.variables[ncvariablename][12:13:,:,:]
fillvalue=ncfile.variables[ncvariablename]._FillValue
ncfile.close()
#combine all data from the chosen variable to make one array for analyses
cdftime=utime('seconds since 1970-01-01 00:00:00')
ncfiletime=cdftime.num2date(TIME)
for i in ncfiletime[:]:
ncfiletime=i
timestr=str(ncfiletime)
d = datetime.strptime(timestr, '%Y-%m-%d %H:%M:%S')
date_string = d.strftime('%H')#combine by same hour
print "the ncfiletime is:", ncfiletime
print "the date_string is:", date_string
if date_string=='00':
time00.append(TSFC00)
elif date_string=='12':
time12.append(TSFC12)
else:
pass
big_arraytime00=N.ma.concatenate(time00)
big_arraytime12=N.ma.concatenate(time12)
MAX00=big_arraytime00.max(axis=0)
MAX12=big_arraytime12.max(axis=0)
#plot output summary stats
map = Basemap(projection='merc',llcrnrlat=-40,urcrnrlat=-33,
llcrnrlon=139.0,urcrnrlon=151.0,lat_ts=0,resolution='i')
map.drawcoastlines()
map.drawstates()
map.readshapefile(shapefile1, 'REGIONS')
x,y=map(*N.meshgrid(LON,LAT))
plottitle='TSFCmax00'
plt.title(plottitle)
CS = map.contourf(x,y,MAX00, ticks, cmap=cmap)
l,b,w,h =0.1,0.1,0.8,0.8
cax = plt.axes([l+w+0.025, b, 0.025, h])
plt.colorbar(CS,cax=cax, drawedges=True)
plt.savefig((os.path.join(OutputFolder, plottitle+'.png')))
plt.show()
plt.close()
map = Basemap(projection='merc',llcrnrlat=-40,urcrnrlat=-33,
llcrnrlon=139.0,urcrnrlon=151.0,lat_ts=0,resolution='i')
map.drawcoastlines()
map.drawstates()
map.readshapefile(shapefile1, 'REGIONS')
x,y=map(*N.meshgrid(LON,LAT))
plottitle='TSFCmax12'
plt.title(plottitle)
CS = map.contourf(x,y,MAX12, ticks, cmap=cmap)
l,b,w,h =0.1,0.1,0.8,0.8
cax = plt.axes([l+w+0.025, b, 0.025, h])
plt.colorbar(CS,cax=cax, drawedges=True)
plt.savefig((os.path.join(OutputFolder, plottitle+'.png')))
plt.show()
plt.close()
_______________________________________________
Tutor maillist - Tutor@python.org
To unsubscribe or change subscription options:
https://mail.python.org/mailman/listinfo/tutor
