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 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文件已传

  • 1
    点赞
  • 9
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值