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)'))