python nc basemap

#!/usr/bin/env python3# -*- coding: utf-8 -*-"""Created on Tue Dec  5 23:10:46 2017@author: chenze"""from mpl_toolkits.basemap import Basemap, cm# requires netcdf4-python (netcdf4-python.googlecode.com)from netCDF4 import Dataset as NetCDFFileimport numpy as npimport matplotlib.pyplot as plt#import h5py# plot rainfall from NWS using special precipitation# colormap used by the NWS, and included in basemap.nc = NetCDFFile('/home/chenze/practice/nws_precip_1day_20171203_netcdf/nws_precip_conus_20061222.nc', mode='r')#nc= NetCDFFile('/home/chenze/practice/nws_precip_1day_20171203_netcdf/nws_precip_1day_20171203_ak.nc')#nc2 = NetCDFFile('/home/chenze/practice/nws_precip_1day_20171203_netcdf/nws_precip_1day_20171203_pr.nc')# data from http://water.weather.gov/precip/print(nc)#,'1\n',nc1,'2\n',nc2)print(nc.variables.keys())prcpvar = nc.variables['amountofprecip']data = 0.01*prcpvar[:]#print(nc.variables['polar_stereographic'][:])print(data)#print(nc2)#print(nc2.variables['percent_of_normal'][:])#print(nc2.variables.keys())#print(nc.variables['normal'][:])#latcorners = nc.variables['lat'][:]#loncorners = -nc.variables['lon'][:]#lon_0 = -nc.variables['true_lon'].getValue()#lat_0 = nc.variables['true_lat'].getValue()# create figure and axes instancesfig = plt.figure(figsize=(8,8))ax = fig.add_axes([0.1,0.1,0.8,0.8])# create polar stereographic Basemap instance.m = Basemap(projection='stere',lon_0=-105,lat_0=90.,lat_ts=60.,\            llcrnrlat=25.,urcrnrlat=50.,\            llcrnrlon=-118.,urcrnrlon=-60.,\            rsphere=6371200.,resolution='l',area_thresh=10000)'''projection  地图投影方式地图底图各参量  lon_0 lat_0 投影中心               llcrnr-- 最小经纬度             urcrnr--最大经纬度  两者对后面计算所有数据的经纬度有影响 其它参量可以看之前推荐的站             '''    # draw coastlines, state and country boundaries, edge of map.m.drawcoastlines()m.drawstates()m.drawcountries()# draw parallels.parallels = np.arange(0.,90,10.)m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)# draw meridiansmeridians = np.arange(180.,360.,10.)m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)ny = data.shape[0]; nx = data.shape[1]print(ny,nx)#经纬度的个数lons, lats = m.makegrid(nx, ny) # get lat/lons of ny by nx evenly space grid.#计算经纬度--根据最大，最小值，以及个数来进行准确计算。print(lons,lats)x, y = m(lons, lats) # compute map proj coordinates.#经纬度坐标，根据地图投影方式，计算对应的系，xy值# draw filled contours.clevs = [0,1,2.5,5,7.5,10,15,20,30,40,50,70,100,150,200,250,300,400,500,600,750]cs = m.contourf(x,y,data,clevs,cmap=cm.s3pcpn)# add colorbar.cbar = m.colorbar(cs,location='bottom',pad="5%")cbar.set_label('mm')# add title#plt.title(prcpvar.long_name+' for period ending '+prcpvar.dateofdata)plt.show()

计算经纬度时，因为用的nc文件改过，所以降水图与地图没有很好的重合。nc文件各数据已给出，懒得改了！
nc文件已传