#!/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 NetCDFFile
import numpy as np
import 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 instances
fig = 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 meridians
meridians = 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文件已传