find location of maximum values
I have a netcdf file that contains hourly temperature data for a whole month. I would like to find the maximum temperature within that file and also the corresponding Latitude and Longitude and Time and then plot this. Below is the code I have so far. I think everything is working except for identifying the corresponding latitude and longitude. The latitude and longitude always seems to be in the same area (from file to file) and they are not where the maximum values are occurring. Any feedback will be greatly appreciated. from netCDF4 import Dataset import matplotlib.pyplot as plt import numpy as N from mpl_toolkits.basemap import Basemap from netcdftime import utime from datetime import datetime from numpy import ma as MA import os TSFCall=[] for (path, dirs, files) in os.walk(MainFolder): for dir in dirs: print dir path=path+'/' for ncfile in files: if ncfile[-3:]=='.nc': print "dealing with ncfiles:", path+ncfile ncfile=os.path.join(path,ncfile) ncfile=Dataset(ncfile, 'r+', 'NETCDF4') TSFC=ncfile.variables['T_SFC'][:] TIME=ncfile.variables['time'][:] LAT=ncfile.variables['latitude'][:] LON=ncfile.variables['longitude'][:] ncfile.close() TSFCall.append(TSFC) big_array=N.ma.concatenate(TSFCall) tmax=MA.max(TSFC) print tmax t=TIME[tmax] indexTIME=N.where(TIME==t) TSFCmax=TSFC[tmax] print t, indexTIME print LAT[tmax], LON[tmax], TSFCmax cdftime=utime('seconds since 1970-01-01 00:00:00') ncfiletime=cdftime.num2date(t) print ncfiletime timestr=str(ncfiletime) d = datetime.strptime(timestr, '%Y-%m-%d %H:%M:%S') date_string = d.strftime('%Y%m%d_%H%M') 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.tissot(LON[tmax], LAT[tmax], 0.15, 100, facecolor='none', edgecolor='black', linewidth=2) map.readshapefile(shapefile1, 'DSE_REGIONS') map.drawcoastlines(linewidth=0.5) map.drawstates() plt.title('test'+' %s'%date_string) ticks=[-5,0,5,10,15,20,25,30,35,40,45,50] CS = map.contourf(x,y,TSFCmax, 15, cmap=plt.cm.jet) l,b,w,h =0.1,0.1,0.8,0.8 cax = plt.axes([l+w+0.025, b, 0.025, h], ) cbar=plt.colorbar(CS, cax=cax, drawedges=True) plt.show() plt.close()
participants (7)
-
"David Köpfer"
-
Aronne Merrelli
-
Benjamin Root
-
Derek Homeier
-
Olivier Delalleau
-
questions anon
-
Vincent Schut