arcgis批量处理nc文件_气象数据处理—

arcgis批量处理nc文件_气象数据处理——nc文件

数据说明

NetCDF(network Common Data Form)网络通用数据格式是一种面向数组型并适于网络共享的数据的描述和编码标准。目前，NetCDF广泛应用于大气科学、水文、海洋学、环境模拟、地球物理等诸多领域。用户可以借助多种方式方便地管理和操作 NetCDF 数据集。

NetCDF全称为network Common Data Format，中文译法为“网络通用数据格式”；netcdf文件开始的目的是用于存储气象科学中的数据，现在已经成为许多数据采集软件的生成文件的格式。

•从数学上来说，netcdf存储的数据就是一个多自变量的单值函数。用公式来说就是f(x,y,z,…)=value；
•函数的自变量x,y,z等在netcdf中叫做维(dimension) 或坐标轴(axix)，
•函数值value在netcdf中叫做变量(Variables).
一个Netcdf文件的结构包括以下对象：
•变量(Variables) ：变量对应着真实的物理数据。
•维(dimension)：一个维对应着函数中的某个自变量，或者说函数图象中的一个坐标轴，在线性代数中就是一个N维向量的一个分量。
•属性(Attribute) ：属性对变量值和维的具体物理含义的注释或者说解释。
（原文链接：CSDN-专业IT技术社区-登录）

Python读取：使用netCDF4的Dataset方法读入文件

# -*- 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网站

odict_keys(['longitude', 'latitude', 'time', 'ssrd', 'ssr', 'fdir', 'strd', 'str'])

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

longitude (1,)
latitude (2,)
time (4392,)
ssrd (4392, 2, 1)
ssr (4392, 2, 1)
fdir (4392, 2, 1)
strd (4392, 2, 1)
str (4392, 2, 1)

#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"))

[1016851 1016852 1016853 1016854 1016855 1016856 1016857 1016858 1016859
1016860]
2016-01-01 19:00:00
2016-01-01 20:00:00
2016-01-01 21:00:00

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

[116.5] [39.875 39.75 ]

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

Python作图：使用matplotlib和Basemap

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个

odict_keys(['longitude', 'latitude', 'time', 't2m'])
(744, 361, 720)

#查看数据经纬度范围，经度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)

0.0 359.5 90.0 -90.0
(720,) (361,)

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

用Basemap画地图

#用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

[[172]] [[72]] [[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)

<class 'numpy.ma.core.MaskedArray'>
(90, 124) (90, 124) (90, 124)

全球温度作图

#全球温度作图，设定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)

<class 'numpy.ma.core.MaskedArray'>
(361, 720) (361, 720) (361, 720)