# -*- coding: utf-8 -*-
from netCDF4 import Dataset
import numpy as np
import sys
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap
from pandas import DataFrame
#数据读入
nc=Dataset('bj2016_18pm.nc')print(nc.variables.keys())#该文件是辐射资料，来自ECMWF网站

#取出各variable的数据看看,数据格式为numpy数组
for var in nc.variables.keys():data=nc.variables[var][:].dataprint(var,data.shape)
#     np.save(var+'.npy',data)

#time variable查看，时间戳变换
#看出是逐时数据
import datetime
time=nc.variables['time'][:].data
print(time[:10])
for i in range(3):tstamp=(time[i]-613608)*3600 #1900年1月1日零时距离1970年1月1日零时有613608个小时date= datetime.datetime.utcfromtimestamp(tstamp)print (date.strftime("%Y-%m-%d %H:%M:%S"))

#查看longitude和 latitude这两个variable
#这个数据比较小，只包含两个格点，所以直接输出了
print(nc.variables['longitude'][:].data,nc.variables['latitude'][:].data)

#查看辐射数据&数据切片
ssrd=nc.variables['ssrd'][:].data
#numpy数组切片，取数
ssrd_timeseq=ssrd[:,1,0]#取出某个格点所有时间的ssrd值
ssrd_timei=ssrd[0,:,:]#取出某个时间点所有格点的ssrd值

nc=Dataset('sfc_201803_2mt.nc')#数据来自ECMWF，为2018年3月份的全球温度的格点数据
print(nc.variables.keys())data=nc.variables['t2m'][:]
print(data.shape)
#d1=data[0,:,:].data##取出某个时刻的温度
#输出结果表明有744个时间序列，经度、纬度的取值有720、361个

#查看数据经纬度范围，经度0-360，其中0~180为西经0~180，180~360为东经0~180；纬度正为北纬，负为南纬
#格点分辨率为0.5度
long= nc.variables['longitude'][:]
lati= nc.variables['latitude'][:]
print(long[0],long[-1],lati[0],lati[-1])
print(long.shape,lati.shape)

#plt对某个时刻的全球温度作图。左半部分为西半球，右边是东半球
#选了某个时间点116作图
plt.contourf(long,lati,data[116,:,:]-273) #转为摄氏度
plt.colorbar()
<matplotlib.colorbar.Colorbar at 0x2b48e65ebe0>

#用Basemap画地图
def graph(lon,lat,target,levelT,colorT,title): b_map=Basemap(resolution='l', area_thresh=10000, projection='cyl', llcrnrlon=min(lon), urcrnrlon=max(lon), llcrnrlat=min(lat),urcrnrlat=max(lat))#llcrnrlon=0, urcrnrlon=360, llcrnrlat=-90,urcrnrlat=90print(type(target))fig=plt.figure(figsize=(9, 6))  #plt.figure(figsize=(12, 8))ax=fig.add_axes([0.1,0.1,0.8,0.8])lon,lat=np.meshgrid(lon,lat)x,y=b_map(lon,lat)print(x.shape,y.shape,target.shape)cs=b_map.contourf(x,y,target,levels=levelT,colors=colorT) #target[0,:,:]  b_map.colorbar(cs)b_map.drawcoastlines(linewidth=1)b_map.drawcountries(linewidth=1.5)plt.title(title,size=20)#plt.savefig('Rainf_0.png',dpi=300)plt.show()plt.close()

#首先查找中国四至范围对应的索引
#4,54为北纬4，54；73+180，133+180为东经73，133
print(np.argwhere(lati==4),np.argwhere(lati==54),
np.argwhere(long==73+180),np.argwhere(long==133+180))
#可以看到，中国范围 对应的纬度索引为72:172,,经度索引为506:626

#对中国范围的温度作图，设定graph的参数
title='2m_temperature'
level_Tair= [-20,-10,0,5,10,15,20,25,30,1000] #[0,2.6,5,8,16,50,100,120,1000]  #[0,2.6,5,8,10,20,25,300,1000]   [0,210,225,240,255,260,300,305,310,1000]
colors = ['#FFFFFF', '#AAF0FF', '#C8DC32',  '#FFBE14', '#FF780A','#FF5A0A', '#F02800',  '#780A00', '#140A00']#注意这里要对经度做变换，原来东半球经度在180~360区间，现在减去180，转为0~180
long=nc.variables['longitude'][506:630]-180#lati= np.flip(nc.variables['latitude'][72:172])
lati= nc.variables['latitude'][62:152]  #[60:164]#datai=np.flipud(data[30,72:172,506:630])-274 #转换为摄氏度
datai=data[30,62:152,506:630]-274 #温度数据切片，选择第30个时间点的温度；将原温度转换为摄氏度graph(long,lati,datai,level_Tair,colors,title)

#全球温度作图，设定graph的参数
title='2m_temperature'
level_Tair= [-20,-10,0,5,10,15,20,25,30,1000] #[0,2.6,5,8,16,50,100,120,1000]  #[0,2.6,5,8,10,20,25,300,1000]   [0,210,225,240,255,260,300,305,310,1000]
colors = ['#FFFFFF', '#AAF0FF', '#C8DC32',  '#FFBE14', '#FF780A','#FF5A0A', '#F02800',  '#780A00', '#140A00']#注意这里要对经度做变换，原来东半球经度在180~360区间，现在减去180，转为0~180
long=nc.variables['longitude'][:]-180#lati= np.flip(nc.variables['latitude'][72:172])
lati= nc.variables['latitude'][:]  #[60:164]#datai=np.flipud(data[30,72:172,506:630])-274 #转换为摄氏度
datai=data[30,:,:]-274 #将原温度转换为摄氏度graph(long,lati,datai,level_Tair,colors,title)