I'm sorry I don't have time to look closely at your code and this may not
be helpful, but just in case... I find it suspicious that you *seem* (by
quickly glancing at the code) to be taking TIME[max(temperature)] instead
of TIME[argmax(temperature)].
-=- Olivier
2011/12/20 questions anon <questions.a...@gmail.com>
> 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()
>
>
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
