import datetime as dt  # Python standard library datetime  moduleimport numpy as npfrom netCDF4 import Dataset  # http://code.google.com/p/netcdf4-python/import matplotlib.pyplot as pltfrom mpl_toolkits.basemap import Basemap, addcyclic, shiftgridnc_f = 'F:/Rpython/lp28/data/air.sig995.2012.nc'  # Your filenamenc_fid = Dataset(nc_f, 'r')  # Dataset is the class behavior to open the file  # and create an instance of the ncCDF4 class# Extract data from NetCDF filelats = nc_fid.variables['lat'][:]  # extract/copy the datalons = nc_fid.variables['lon'][:]time = nc_fid.variables['time'][:]air = nc_fid.variables['air'][:]  # shape is time, lat, lon as shown abovetime_idx = 237  # some random day in 2012# Python and the renalaysis are slightly off in time so this fixes that problemoffset = dt.timedelta(hours=48)# List of all times in the file as datetime objectsdt_time = [dt.date(1, 1, 1) + dt.timedelta(hours=t) - offset\           for t in time]cur_time = dt_time[time_idx]# Plot of global temperature on our random day#fig = plt.figure()#fig.subplots_adjust(left=0., right=1., bottom=0., top=0.9)fig=plt.figure(figsize=(16,9))ax=fig.add_subplot(111)# Setup the map. See http://matplotlib.org/basemap/users/mapsetup.html# for other projections.m=Basemap(llcrnrlat=-90,urcrnrlat=90,llcrnrlon=-180,urcrnrlon=180,resolution='c')m.drawparallels(np.arange(-90, 90 + 1, 30), labels = [1, 0, 0, 0],fontsize=14,linewidth='0.2',color='black')m.drawmeridians(np.arange(-180, 180 + 1, 60), labels = [0, 0, 0, 1],fontsize=14,linewidth='0.2',color='black')m.drawcoastlines()m.drawmapboundary()air2=air[time_idx,:,:]# 将0~360转化为-180~180# Make the plot continuousair_cyclic, lons_cyclic = addcyclic(air2, lons)# Shift the grid so lons go from -180 to 180 instead of 0 to 360.air_cyclic, lons_cyclic = shiftgrid(180., air_cyclic, lons_cyclic, start=False) #False True# Create 2D lat/lon arrays for Basemaplon2d, lat2d = np.meshgrid(lons_cyclic, lats)# Transforms lat/lon into plotting coordinates for projectionx, y = m(lon2d, lat2d)# 查看数组纬度dimen = np.array(x).shapeprint(dimen)dimen1 = np.array(y).shapeprint(dimen1)dimen2 = np.array(air_cyclic).shapeprint(dimen2)# Plot of air temperature with 11 contour intervalsclevs = np.arange(200,320,10)cs = m.contourf(x,y,air_cyclic,clevs,cmap='gist_rainbow',extend='both')cbar=plt.colorbar(cs,shrink=0.75,orientation='vertical',extend='both',pad=0.025,aspect=20) #orientation='horizontal'#cbar = plt.colorbar(cs, orientation='horizontal', shrink=0.5)cbar.set_label("%s (%s)" % (nc_fid.variables['air'].var_desc,nc_fid.variables['air'].units))plt.title("%s on %s" % (nc_fid.variables['air'].var_desc, cur_time))plt.savefig('F:/Rpython/lp28/plot33.png',dpi=800)plt.show()

darwin = {'name': 'Darwin, Australia', 'lat': -12.45, 'lon': 130.83}# Find the nearest latitude and longitude for Darwinlat_idx = np.abs(lats - darwin['lat']).argmin()lon_idx = np.abs(lons - darwin['lon']).argmin()# Simple example: temperature profile for the entire year at Darwin.# Open a new NetCDF file to write the data to. For format, you can choose from# 'NETCDF3_CLASSIC', 'NETCDF3_64BIT', 'NETCDF4_CLASSIC', and 'NETCDF4'w_nc_fid = Dataset('F:/Rpython/lp28/data/darwin_2012.nc', 'w', format='NETCDF4')w_nc_fid.description = "NCEP/NCAR Reanalysis %s from its value at %s. %s" %\                      (nc_fid.variables['air'].var_desc.lower(),darwin['name'], nc_fid.description)# Using our previous dimension info, we can create the new time dimension# Even though we know the size, we are going to set the size to unknownw_nc_fid.createDimension('time', None)w_nc_dim = w_nc_fid.createVariable('time', nc_fid.variables['time'].dtype,('time',))# You can do this step yourself but someone else did the work for us.for ncattr in nc_fid.variables['time'].ncattrs():    w_nc_dim.setncattr(ncattr, nc_fid.variables['time'].getncattr(ncattr))# Assign the dimension data to the new NetCDF file.w_nc_fid.variables['time'][:] = timew_nc_var = w_nc_fid.createVariable('air', 'f8', ('time'))w_nc_var.setncatts({'long_name': u"mean Daily Air temperature",\                    'units': u"degK", 'level_desc': u'Surface',\                    'var_desc': u"Air temperature",\                    'statistic': u'Mean\nM'})w_nc_fid.variables['air'][:] = air[time_idx, lat_idx, lon_idx]w_nc_fid.close()  # close the new file# A plot of the temperature profile for Darwin in 2012fig = plt.figure()plt.plot(dt_time, air[:, lat_idx, lon_idx], c='r')plt.plot(dt_time[time_idx], air[time_idx, lat_idx, lon_idx], c='b', marker='o')plt.text(dt_time[time_idx], air[time_idx, lat_idx, lon_idx], cur_time,ha='right')fig.autofmt_xdate()plt.ylabel("%s (%s)" % (nc_fid.variables['air'].var_desc,nc_fid.variables['air'].units))plt.xlabel("Time")plt.title("%s from\n%s for %s" % (nc_fid.variables['air'].var_desc,darwin['name'], cur_time.year))plt.savefig('F:/Rpython/lp28/plot33.2.png',dpi=800)plt.show()