MeteoInfoLab脚本示例:创建netCDF文件(合并文件)

在MeteoInfoLab中增加了创建netCDF文件并写入数据的功能,这里利用合并多个netCDF文件为一个新的netCDF文件为例。

1、创建一个可写入的netCDF文件对象(下面用ncfile表示),用addfile函数,第一个参数是文件名,第二次参数‘c‘表示创建新的netCDF文件。
ncfile = addfile(outfn, ‘c‘)

2、添加维(Dimensions),用ncfile的adddim函数,两个参数分别是维名称和维长度。
stn = 26564
stdim = ncfile.adddim(‘station‘, stn)

3、添加全局属性,用ncfile的addgroupattr函数,两个参数分别是属性名称和属性值。
ncfile.addgroupattr(‘Conventions‘, ‘Unidata Observation Dataset v1.0‘)

4、添加变量,用ncfile的addvar函数,三个参数分别是变量名称、变量数据类型和变量维列表。并给变量添加一些属性用变量对象的addattr函数。
var = ncfile.addvar(‘streamflow‘, ‘float‘, [tdim, stdim])
var.addattr(‘long_name‘, ‘River Flow‘)
var.addattr(‘units‘, ‘meter^3 / sec‘)

5、创建netCDF文件,上面的步骤都是对netCDF文件的设置,完成之后需要create函数表示设置已完成,可以开始写数据了。
ncfile.create()

6、写数据,用ncfile的write函数,参数分别是变量和要写入的数组,origin是各维起点的列表,可以控制写入某个变量的部分数据。如果要写入的数组的维数和变量的维数不一致,需要用reshape函数来确保二者有相同的维数。
ncfile.write(variables[3], flow, origin=origin)

7、关闭netCDF文件,数据写完之后需要关闭netCDF文件。
ncfile.flush()
ncfile.close()

示例脚本程序(合并4个不同时次的netCDF文件为一个新的netCDF文件):

datadir = ‘D:/Temp/nc‘
outfn = os.path.join(datadir, ‘join_file.nc‘)
#New netCDF file
ncfile = addfile(outfn, ‘c‘)
#Add dimensions
stn = 26564
recdim = ncfile.adddim(‘recNum‘, stn)
stdim = ncfile.adddim(‘station‘, stn)
iddim = ncfile.adddim(‘id_len‘, 11)
tdim = ncfile.adddim(‘time‘, 4)
#Add global attributes
ncfile.addgroupattr(‘Conventions‘, ‘Unidata Observation Dataset v1.0‘)
ncfile.addgroupattr(‘cdm_datatype‘, ‘Station‘)
ncfile.addgroupattr(‘geospatial_lat_max‘, ‘90.0‘)
ncfile.addgroupattr(‘geospatial_lat_min‘, ‘-90.0‘)
ncfile.addgroupattr(‘geospatial_lon_max‘, ‘180.0‘)
ncfile.addgroupattr(‘geospatial_lon_min‘, ‘-180.0‘)
ncfile.addgroupattr(‘stationDimension‘, ‘station‘)
ncfile.addgroupattr(‘missing_value‘, -8.9999998E15)
ncfile.addgroupattr(‘stream_order_output‘, 1)
#Add variables
variables = []
var = ncfile.addvar(‘latitude‘, ‘float‘, [stdim])    #Latitude
var.addattr(‘long_name‘, ‘station latitude‘)
var.addattr(‘units‘, ‘degrees_north‘)
variables.append(var)
var = ncfile.addvar(‘longitude‘, ‘float‘, [stdim])    #Longitude
var.addattr(‘long_name‘, ‘station longitude‘)
var.addattr(‘units‘, ‘degrees_east‘)
variables.append(var)
var = ncfile.addvar(‘altitude‘, ‘float‘, [stdim])    #Altitude
var.addattr(‘long_name‘, ‘station altitude‘)
var.addattr(‘units‘, ‘meters‘)
variables.append(var)
var = ncfile.addvar(‘streamflow‘, ‘float‘, [tdim, stdim])    #Stream flow - Add time dimension
var.addattr(‘long_name‘, ‘River Flow‘)
var.addattr(‘units‘, ‘meter^3 / sec‘)
variables.append(var)
tvar = ncfile.addvar(‘time‘, ‘int‘, [tdim])
tvar.addattr(‘long_name‘, ‘time‘)
tvar.addattr(‘units‘, ‘hours since 1900-01-01 00:00:0.0‘)
#Creat netCDF file
ncfile.create()
#Write data
stime = datetime.datetime(2015,8,2,0)
etime = datetime.datetime(2015,8,2,3)
st = datetime.datetime(1900,1,1)
fi = 0
while stime <= etime:
    print stime
    fn = os.path.join(datadir, stime.strftime(‘%Y%m%d%H‘) + ‘00.CHRTOUT_DOMAIN2‘)
    if os.path.exists(fn):
        print ‘\t‘ + fn
        f = addfile(fn)
        hours = (stime - st).total_seconds() // 3600
        origin = [fi]
        ncfile.write(tvar, array([hours]), origin=origin)
        if fi == 0:
            lat = f[‘latitude‘][:]
            ncfile.write(variables[0], lat)
            lon = f[‘longitude‘][:]
            ncfile.write(variables[1], lon)
            alt = f[‘altitude‘][:]
            ncfile.write(variables[2], alt)
        flow = f[‘streamflow‘][:]
        origin = [fi, 0]
        shape = [1, stn]
        flow = flow.array.reshape(shape)
        ncfile.write(variables[3], flow, origin=origin)
        fi += 1
    stime = stime + datetime.timedelta(hours=1)

#close netCDF file
ncfile.flush()
ncfile.close()

print ‘Finished!‘

对合并后的文件绘图:

f = addfile(‘D:/Temp/nc/join_file.nc‘)
lon = f[‘longitude‘][:]
lat = f[‘latitude‘][:]
var = f[‘streamflow‘]
flow = var[1,:]
axesm()
mlayer = shaperead(‘D:/Temp/Map/bou2_4p.shp‘)
geoshow(mlayer)
levs = arange(0, 0.1, 0.005)
layer = scatterm(lon, lat, flow, levs, edge=False)
colorbar(layer)
t = f.gettime(1)
title(‘River Flow (‘ + t.strftime(‘%Y-%m-%d %Hh)‘))

时间: 2024-12-26 15:44:56

MeteoInfoLab脚本示例:创建netCDF文件(合并文件)的相关文章

MeteoInfoLab脚本示例:加载图片和透明图层

MeteoInfoLab的georead函数提供了读取shape文件.image文件(JPG.PNG等,需要有相应的地理定位文件)文件生成图层的功能(事实上shaperead也是同样的功能,不过函数名太过具体).对于Polygon图层可以设置其颜色透明度,可以用makecolors函数生成颜色列表,函数中设置参数alpha(255是不透明,0是完全透明).在contourfm函数中设置edgecolor参数可以绘制颜色区域的边框. 脚本程序: f = addfile('D:/Temp/GrADS

MeteoInfoLab脚本示例:MODIS AOD

MODIS的气溶胶光学厚度(AOD)产品应用很广,数据可以在Giovanni上下载:http://disc.sci.gsfc.nasa.gov/giovanni/overview/index.html.有HDF, netCDF和文本格式,这里示例打开netCDF格式数据,提取AOD数据并绘图. 脚本程序如下: f = addfile('D:/Temp/nc/MOD08_D3.A2015121.051.2015122103938.pscs_000500931513.Optical_Depth_La

MeteoInfoLab脚本示例:计算垂直螺旋度

尝试编写MeteoInfoLab脚本计算垂直螺旋度,结果未经验证. 脚本程序: print 'Open data files...' f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc') f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc') f_omega = addfile('D:/Temp/nc/omega.2011.nc') print 'Calculate vertical helicity...' tidx = 173

MeteoInfoLab脚本示例:FY-3C全球火点HDF数据

FY-3C全球火点HDF数据包含一个FIRES二维变量,第一维是火点数,第二维是一些属性,其中第3.4列分别是火点的纬度和经度.下面的脚本示例读出所有火点经纬度并绘图. 脚本程序: #Add data file fn = 'D:/Temp/hdf/FY3C_VIRRX_GBAL_L2_GFR_MLT_GLL_20150811_POAD_1000M_MS.HDF' f = addfile(fn) #Get data variable v = f['FIRES'] #Get data array d

MeteoInfoLab脚本示例:合并数组

对于全球数据来说,经度要么是-180 - 180,要么是0 - 360,都会存在边界数据不连续的问题.比如0 - 360的数据,怎么得到 -20 - 30度的连续格点数据就是个问题(跨越了数据的经度边界),在MeteoInfoLab中可以用DimArray或MIArray的join函数来将两个数组合并为一个,参数分别是另一个数组和合并的维的序号,比如下面例子中的二维数组(y, x),经度维是第二维(序号为1),也就是说两个数组沿着经度维合并.还需要给合并后的数组的经度维设置新的经度值. 脚本程序

MeteoInfoLab脚本示例:加载地图图层

应用最广泛的的地图数据应该是shape格式,网络上有很多免费下载资源.MeteoInfoLab中读取shape文件的函数是shaperead,参数即文件名,返回数据包含图形和属性信息的图层对象.矢量图层按照图元类型通常分点.线.面三种,都可以通过geoshow函数加入到地图坐标系(Axesm)中.在geoshow函数中可以指定图层的显示的Symbol(颜色等),可以设置facecolor, edgecolor, size等,这种情况下图层的所有图元都会以同一个Symbol来显示.也可以利用图层的

MeteoInfoLab脚本示例:数据投影-FLEXPART

FLEXPART是一个类似HYSPLIT的扩散模式,它输出的netcdf文件参照了WRF,可惜全局属性没有写全,比如只有一个投影名称(例如Lambert),没有相关的投影参数:中央经度,标准纬度等等.必须查阅WRF的头文件才能重建投影(为什么不照猫画虎把属性写全呢?).数据的经纬度坐标是有的,但在Lambert投影下的坐标没有,可以通过projectxy函数获得投影下的x, y坐标,其中的lon, lat是数据左下角的经纬度. 脚本程序: f = addfile('D:/Temp/nc/head

MeteoInfoLab脚本示例:FY-3A AOD HDF数据

FY3A卫星有HDF格式的AOD产品数据,全球范围,分辨率为0.05度.读取数据文件变量后要重新设定X/Y维,数据是Y轴反向的,且需要除以10000得到AOD值. 脚本程序: #Add data file fn = 'D:/Temp/hdf/FY3A_MERSI_GBAL_L3_ASL_MLT_GLL_20140331_AOAM_5000M_MS.HDF' f = addfile(fn) #Get data variable v = f['Aerosol_Optical_Thickness_of

MeteoInfoLab脚本示例:获取气团轨迹每个节点的气象数据

读取HYSPLIT输出的轨迹数据文件和相应时间的气象数据文件,生成轨迹图层,循环每条轨迹的节点,读出该节点的经度.纬度.气压.时间,通过对气象数据插值获得该节点的气象数据. 脚本程序: #----------------------------------------------------- # Author: Yaqiang Wang # Date: 2015-9-30 # Purpose: Get meteorological data along trajectory # Note: S