上一篇介绍了配置PALM运行的基础，这一篇介绍创建静态驱动模型来自定义的地形。
几点说明：
1）输入的文件格式为netcdf
2）采用python编写
地形
这是我重采样后的地形图片，200像素*100像素。制作步骤如下：
1）下载相应位置的影像
2）要素分类
3）重采样
地形编写netcdf
代码声明
1）基础代码下载自PALM官方网站，其余为笔者编写（转载请声明谢谢）
2）我采用的方法是根据像素值来确定位置（我一直在寻找方法怎么样快速创建netcdf，目前没有找到比较好的方法，如果有好方法请不吝赐教）
3）代码编写需要参考官网：https://palm.muk.uni-hannover.de/trac/wiki/doc/app/iofiles/pids/static
#------------------------------------------------------------------------------#
# This file is part of the PALM model system.
#
# PALM is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# PALM is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with
# PALM. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright 1997-2019 Leibniz Universitaet Hannover
#------------------------------------------------------------------------------#
import math
from netCDF4 import Dataset
import numpy as np
import datetime
class StaticDriver:
""" This is an example script to generate static drivers for PALM.
You can use it as a starting point for creating your setup specific
driver.
"""
def __init__(self):
""" Open the static driver as NetCDF4 file. Here, you have to give the
full path to the static driver that shall be created. Existing file
with same name is deleted.
"""
print('Opening file...')
self.nc_file = Dataset('ysqPalm_static.nc', 'w', format='NETCDF4')
def write_global_attributes(self):
""" Write global attributes to static driver. """
print("Writing global attributes...")
# mandatory global attributes
self.nc_file.origin_lon = 118.7852552 # used to initialize coriolis parameter
self.nc_file.origin_lat = 31.9140447 # (overwrite initialization_parameters)
self.nc_file.origin_time = '2018-07-19 07:00:00 +08'
self.nc_file.origin_x = 0
self.nc_file.origin_y = 0
self.nc_file.origin_z = 0
self.nc_file.rotation_angle = 0.0
# optional global attributes
self.nc_file.creation_date = str(datetime.datetime.now())
def define_dimensions(self): #定义维度
""" Set dimensions on which variables are defined. """
print("Writing dimensions...")
# specify general grid parameters
# these values must equal those set in the initialization_parameters
self.nx = 199
self.ny = 99
self.nz = 30
dx = 1.0
dy = 1.0
dz = 1.0
# create soil grid (only relevant if land surface module is used)
dz_soil = np.array((0.01, 0.02, 0.04, 0.06, 0.14, 0.26, 0.54, 1.86))
zsoil_fullLayers = np.zeros_like(dz_soil)
zsoil_fullLayers = np.around([np.sum(dz_soil[:zs]) for zs in np.arange(1,len(dz_soil)+1)],2)
zsoil_array = zsoil_fullLayers - dz_soil/2.
# create plant canopy grid (only relevant of plant canopy module is used)
zlad_array = np.arange(0,10,1)
# mandatory dimensions
self.nc_file.createDimension('x' ,self.nx+1)
self.x = self.nc_file.createVariable('x', 'f8', ('x',))
self.x.units = 'm'
self.x.standard_name = 'x coordinate of cell centers'
self.x[:] = np.arange(0,(self.nx+1)*dx,dx)
# self.nc_file.createDimension('xu', self.nx)
# self.xu = self.nc_file.createVariable('xu', 'f8', ('xu',))
# self.xu.units = 'm'
# self.xu.standard_name = 'x coordinate of cell edges'
# self.xu[:] = np.arange(dx/2.,self.nx*dx,dx)
self.nc_file.createDimension('y', self.ny+1)
self.y = self.nc_file.createVariable('y', 'f8', ('y',))
self.y.units = 'm'
self.y.standard_name = 'y coordinate of cell centers'
self.y[:] = np.arange(0,(self.ny+1)*dy,dy)
# self.nc_file.createDimension('yv', self.ny)
# self.yv = self.nc_file.createVariable('yv', 'f8', ('yv',))
# self.yv.units = 'm'
# self.yv.standard_name = 'y coordinate of cell edges'
# self.yv[:] = np.arange(dy/2.,self.ny*dy,dy)
# if your simulation uses a stretched vertical grid, you need to
# modify the z and zw coordinates,
# e.g. z_array = (...) and zw_array = (...)
z_array = np.array(
[0,1,3,5,7,9,11,13,15]
)
self.nc_file.createDimension('z',9)
self.z = self.nc_file.createVariable('z', 'f8', ('z',))
self.z.units = 'm'
self.z.standard_name = 'z coordinate of cell centers'
self.z[:] = z_array
# zw_array = np.arange(0,(self.nz+2)*dz,dz)
# self.nc_file.createDimension('zw', self.nz+2)
# self.zw = self.nc_file.createVariable('zw', 'f8', ('zw',))
# self.zw.units = 'm'
# self.zw.standard_name = 'z coordinate of cell edges'
# self.zw[:] = zw_array
# optional dimensions, uncomment if needed
self.nc_file.createDimension('zsoil', len(zsoil_array))
self.zsoil = self.nc_file.createVariable('zsoil', 'f8', ('zsoil',))
self.zsoil.positive = 'down'
self.zsoil.units = 'm'
self.zsoil.standard_name = 'depth below land'
self.zsoil[:] = zsoil_array
#print(self.zsoil)
#self.nc_file.createDimension('zlad', len(zlad_array))
#self.zlad = self.nc_file.createVariable('zlad', 'f8', ('zlad',))
#self.zlad.units = 'm'
#self.zlad.standard_name = 'z coordinate of resolved plant canopy'
#self.zlad[:] = zlad_array
#self.nc_file.createDimension('nalbedo_pars', 3)
#self.nalbedo_pars = self.nc_file.createVariable(
#'nalbedo_pars', 'i1', ('nalbedo_pars',))
#self.nalbedo_pars[:] = np.arange(0,3)
# self.nc_file.createDimension('nbuilding_pars', 46)
# self.nbuilding_pars = self.nc_file.createVariable(
# 'nbuilding_pars', 'i4', ('nbuilding_pars',))
# self.nbuilding_pars[:] = np.arange(0,46)
#self.nc_file.createDimension('nsoil_pars', 8)
#self.nsoil_pars = self.nc_file.createVariable(
#'nsoil_pars', 'i1', ('nsoil_pars',))
#self.nsoil_pars[:] = np.arange(0,8)
# self.nc_file.createDimension('nsurface_fraction', 3)
# self.nsurface_fraction = self.nc_file.createVariable(
# 'nsurface_fraction', 'i1', ('nsurface_fraction',))
# self.nsurface_fraction[:] = np.arange(0,3)
#self.nc_file.createDimension('nvegetation_pars', 12)
#self.nvegetation_pars = self.nc_file.createVariable(
#'nvegetation_pars', 'i1', ('nvegetation_pars',))
#self.nvegetation_pars[:] = np.arange(0,12)
self.nc_file.createDimension('npavement_pars', 4)
self.npavement_pars = self.nc_file.createVariable(
'npavement_pars', 'i4', ('npavement_pars',))
self.npavement_pars[:] = np.arange(0,4)
self.nc_file.createDimension('npavement_subsurface_pars', 2)
self.npavement_subsurface_pars = self.nc_file.createVariable(
'npavement_subsurface_pars', 'i4', ('npavement_subsurface_pars',))
self.npavement_subsurface_pars[:] = np.arange(0,2)
self.nc_file.createDimension('nwater_pars', 7)
self.nwater_pars = self.nc_file.createVariable(
'nwater_pars', 'i1', ('nwater_pars',))
self.nwater_pars[:] = np.arange(0,7)
def add_variables(self):
""" Uncomment variables below as you like.
Be aware that some variables depend on others. For a description of
each variable, please have a look at the documentation at
palm.muk.uni-hannover.de/trac/wiki/doc/app/iofiles/pids/static
An example of how you modify the variables is given below:
building_2d_array = np.ones((self.ny+1,self.nx+1)) * -9999.9
south_wall, north_wall, left_wall, right_wall = 20, 25, 20, 25
building_2d_array[south_wall:north_wall,left_wall:right_wall] = 50
nc_buildings_2d = self.nc_file.createVariable(
'buildings_2d', 'f4', ('y','x'),fill_value=-9999.9)
nc_buildings_2d.lod = 1
nc_buildings_2d[:,:] = building_2d_array
"""
print("Writing variables...")
array_image_txt = np.loadtxt('PALM_AREA.bmp.txt')
row_array_image = array_image_txt.shape[0]
col_array_image = array_image_txt.shape[1]
print(row_array_image,col_array_image)
array_image = np.ones((100,200))
for i in range(row_array_image):
for j in range(col_array_image):
array_image[i,j] = array_image_txt[99-i,j]
## topography variables
building_2d_array = np.ones((self.ny + 1, self.nx + 1)) * -9999.9
for i in range(row_array_image):
for j in range(col_array_image):
if array_image[i,j] == 228.0:
building_2d_array[i,j] = 16.0
if array_image[i,j] == 196.0:
building_2d_array[i,j] = 8.0
nc_buildings_2d = self.nc_file.createVariable(
'buildings_2d', 'f4', ('y','x'), fill_value=-9999.9)
nc_buildings_2d.long_name = "buildings"
nc_buildings_2d.res_orig = 1.0
nc_buildings_2d.units = "m"
#nc_buildings_2d.FillValue = -9999.9
nc_buildings_2d.coordinates = "E_UTM N_UTM lon lat"
nc_buildings_2d.lod = 1
nc_buildings_2d[:,:] = building_2d_array
#nc_buildings_3d = self.nc_file.createVariable(
#'buildings_3d', 'i1', ('z','y','x'),fill_value=-127)
#nc_buildings_3d.lod = 2
#nc_buildings_3d[:,:,:] = -127
building_id_array = np.ones((self.ny + 1, self.nx + 1)) * -9999
for i in range(row_array_image):
for j in range(col_array_image):
if array_image[i, j] == 228.0:
building_id_array[i, j] = 1
if array_image[i, j] == 196.0:
building_id_array[i, j] = 2
nc_building_id = self.nc_file.createVariable(
'building_id', 'i4', ('y','x'), fill_value=-9999)
nc_building_id.long_name = "building_id"
nc_building_id.res_orig = 1.0
#nc_building_id._FillValue = -9999
nc_building_id.coordinates = "E_UTM N_UTM lon lat"
nc_building_id[:,:] = building_id_array
zt_array = np.ones((self.ny + 1, self.nx + 1)) * -9999.9
for i in range(row_array_image):
for j in range(col_array_image):
zt_array[i,j] = 0.0
nc_zt = self.nc_file.createVariable(
'zt', 'f4', ('y','x'), fill_value=-9999.9)
nc_zt.long_name = 'orography'
nc_zt.res_orig = 1.0
nc_zt.units = 'm'
#nc_zt._FillValue = -9999.9
nc_zt.coordinate = "E_UTM N_UTM lon lat"
nc_zt[:,:] = zt_array
## surface variables
#nc_albedo_pars = self.nc_file.createVariable(
#'albedo_pars', 'f4', ('nalbedo_pars','y','x'), fill_value=-9999.9)
#nc_albedo_pars.long_name = 'albedo parameters'
#nc_albedo_pars[:,:,:] = -9999.9
#nc_albedo_type = self.nc_file.createVariable(
#'albedo_type', 'f4', ('y','x'), fill_value=-9999.9)
#nc_albedo_type[:,:] = -9999.9
# nc_building_pars = self.nc_file.createVariable(
# 'building_pars', 'f4', ('nbuilding_pars','y','x'), fill_value=-9999.9)
# nc_building_pars.long_name = 'building parameters'
# nc_building_pars[:,:,:] = -9999.9
building_type_array = np.ones((self.ny + 1, self.nx + 1)) * -127
for i in range(row_array_image):
for j in range(col_array_image):
if array_image[i, j] == 228.0:
building_type_array[i, j] = 3
if array_image[i, j] == 196.0:
building_type_array[i, j] = 3
nc_building_type = self.nc_file.createVariable(
'building_type', 'i1', ('y','x'), fill_value=-127)
nc_building_type.long_name = "building type"
nc_building_type.res_orig = 1.0
#nc_building_type._FillValue = -127
nc_building_type.coordinates = "E_UTM N_UTM lon lat"
nc_building_type[:,:] = building_type_array
# nc_pavement_pars = self.nc_file.createVariable(
# 'pavement_pars', 'f4', ('npavement_pars','y','x'), fill_value=-9999.9)
# nc_pavement_pars_array0 = np.ones((self.ny + 1, self.nx + 1)) * -9999.9
# for i in range(row_array_image):
# for j in range(col_array_image):
# if array_image[i, j] == 0:
# nc_pavement_pars_array0[i,j] = 1.0*math.pow(10,-2)
# if array_image[i,j] == 78:
# nc_pavement_pars_array0[i,j] = 1.0*math.pow(10,-2)
# if array_image[i,j] == 158:
# nc_pavement_pars_array0[i,j] = 1.0*math.pow(10,-2)
# nc_pavement_pars_array1 = np.ones((self.ny + 1, self.nx + 1)) * -9999.9
# for i in range(row_array_image):
# for j in range(col_array_image):
# if array_image[i, j] == 0:
# nc_pavement_pars_array1[i,j] = 1.0*math.pow(10,-4)
# if array_image[i,j] == 78:
# nc_pavement_pars_array1[i,j] = 1.0*math.pow(10,-4)
# if array_image[i,j] == 158:
# nc_pavement_pars_array1[i,j] = 1.0*math.pow(10,-4)
# nc_pavement_pars_array2 = np.ones((self.ny + 1, self.nx + 1)) * -9999.9
# for i in range(row_array_image):
# for j in range(col_array_image):
# if array_image[i, j] == 0:
# nc_pavement_pars_array2[i,j] = 18
# if array_image[i,j] == 78:
# nc_pavement_pars_array2[i,j] = 20
# if array_image[i,j] == 158:
# nc_pavement_pars_array2[i,j] = 11
# nc_pavement_pars_array3 = np.ones((self.ny + 1, self.nx + 1)) * -9999.9
# for i in range(row_array_image):
# for j in range(col_array_image):
# if array_image[i, j] == 0:
# nc_pavement_pars_array3[i,j] = 0.95
# if array_image[i,j] == 78:
# nc_pavement_pars_array3[i,j] = 0.95
# if array_image[i,j] == 158:
# nc_pavement_pars_array3[i,j] = 0.95
#
# nc_pavement_pars[0,:,:] = nc_pavement_pars_array0
# nc_pavement_pars[1, :, :] = nc_pavement_pars_array1
# nc_pavement_pars[2, :, :] = nc_pavement_pars_array2
# nc_pavement_pars[3, :, :] = nc_pavement_pars_array3
# #nc_pavement_pars[:,:,:] = -9999.9
#
# nc_pavement_subsurface_pars = self.nc_file.createVariable(
# 'pavement_subsurface_pars', 'f4',
# ('npavement_subsurface_pars','zsoil','y','x'), fill_value=-9999.9)
# nc_pavement_subsurface_pars_array0 = np.ones((8,self.ny + 1, self.nx + 1)) * -9999.9
# for i in range(row_array_image):
# for j in range(col_array_image):
# if array_image[i, j] == 0:
# nc_pavement_subsurface_pars_array0[:,i,j] = 1
# if array_image[i,j] == 78:
# nc_pavement_subsurface_pars_array0[:,i,j] = 1
# if array_image[i,j] == 158:
# nc_pavement_subsurface_pars_array0[:,i,j] = 1
# nc_pavement_subsurface_pars_array1 = np.ones((8,self.ny + 1, self.nx + 1)) * -9999.9
# for i in range(row_array_image):
# for j in range(col_array_image):
# if array_image[i, j] == 0: #1.94*math.pow(10,6)
# nc_pavement_subsurface_pars_array1[:,i,j] = 1.94*math.pow(10,6)
# if array_image[i,j] == 78:
# nc_pavement_subsurface_pars_array1[:,i,j] = 1.94*math.pow(10,6)
# if array_image[i,j] == 158:
# nc_pavement_subsurface_pars_array1[:,i,j] = 1.94*math.pow(10,6)
# nc_pavement_subsurface_pars[0,:,:,:] = nc_pavement_subsurface_pars_array0
# nc_pavement_subsurface_pars[1, :, :, :] = nc_pavement_subsurface_pars_array1
# #nc_pavement_subsurface_pars[:,:,:,:] = -9999.9
pavement_type_array = np.ones((self.ny + 1, self.nx + 1)) * -127
for i in range(row_array_image):
for j in range(col_array_image):
if array_image[i, j] == 0:
pavement_type_array[i,j] = 1
if array_image[i,j] == 78:
pavement_type_array[i,j] = 2
if array_image[i,j] == 158:
pavement_type_array[i,j] = 3
nc_pavement_type = self.nc_file.createVariable(
'pavement_type', 'i1', ('y','x'), fill_value=-127)
nc_pavement_type.long_name = "pavement type"
nc_pavement_type.res_orig = 1.0
#nc_pavement_type._FillValue = -127
nc_pavement_type.coordinates = "E_UTM N_UTM lon lat"
nc_pavement_type[:,:] = pavement_type_array
#nc_soil_pars = self.nc_file.createVariable(
#'soil_pars', 'f', ('nsoil_pars','zsoil','y','x'), fill_value=-9999.9)
#nc_soil_pars.long_name = 'soil parameters'
#nc_soil_pars.lod = 2 # if set to 1, adjust dimensions of soil_pars
#nc_soil_pars[:,:,:,:] = -9999.9
soil_type_array = np.ones((self.ny + 1, self.nx + 1)) * -127
for i in range(row_array_image):
for j in range(col_array_image):
soil_type_array[i,j] = 2
nc_soil_type = self.nc_file.createVariable(
'soil_type', 'i1', ('y','x'), fill_value=-127)
nc_soil_type.long_name = "soil type"
nc_soil_type.res_orig = 1.0
#nc_soil_type._FillValue = -127
nc_building_type.coordinates = "E_UTM N_UTM lon lat"
nc_soil_type.lod = 1 # if set to 1, adjust dimensions of soil_type
nc_soil_type[:,:] = soil_type_array
## required, if several surface types are defined at the same location
# surface_fraction_array = np.ones((self.ny + 1, self.nx + 1)) * -9999.9
# for i in range(row_array_image):
# for j in range(col_array_image):
# if array_image[i, j] == 0.447:
# surface_fraction_array[i,j] = 1
# nc_surface_fraction = self.nc_file.createVariable(
# 'surface_fraction', 'f', ('nsurface_fraction','y','x'), fill_value=-9999.9)
# nc_surface_fraction[0,:,:] = surface_fraction_array # vegetation fraction
# nc_surface_fraction[1,:,:] = 0 # pavement fraction
# nc_surface_fraction[2,:,:] = 0 # water fraction
#nc_vegetation_pars = self.nc_file.createVariable(
#'vegetation_pars', 'f', ('nvegetation_pars','y','x'), fill_value=-9999.9)
#nc_vegetation_pars.long_name = 'vegetation parameters'
#nc_vegetation_pars[:,:,:] = -9999.9
vegetation_type_array = np.ones((self.ny + 1, self.nx + 1)) * -127
for i in range(row_array_image):
for j in range(col_array_image):
if array_image[i,j] == 112.0:
vegetation_type_array[i,j] = 6
if array_image[i,j] == 15.0:
vegetation_type_array[i,j] = 6
if array_image[i,j] == 24.0:
vegetation_type_array[i,j] = 15
if array_image[i,j] == 120.0:
vegetation_type_array[i,j] = 3
if array_image[i,j] == 248.0:
vegetation_type_array[i,j] = 1
nc_vegetation_type = self.nc_file.createVariable(
'vegetation_type', 'i1', ('y','x'), fill_value=-127)
nc_vegetation_type.long_name = "vegetation type"
nc_vegetation_type.res_orig = 1.0
#nc_vegetation_type._FillValue = -127
nc_vegetation_type.coordinates = "E_UTM N_UTM lon lat"
nc_vegetation_type[:,:] = vegetation_type_array
#nc_street_type = self.nc_file.createVariable(
#'street_type', 'i1', ('y','x'),
#fill_value=-127)
#nc_street_type[:,:] = -127
water_par = np.ones((self.ny + 1, self.nx + 1)) * -9999.9
water_pars_array = np.ones((self.ny + 1, self.nx + 1)) * -9999.9
nc_water_pars = self.nc_file.createVariable(
'water_pars', 'f4', ('nwater_pars','y','x'), fill_value=-9999.9)
for i in range(row_array_image):
for j in range(col_array_image):
if array_image[i,j] == 16.0:
water_pars_array[i,j] = 298.35
for i in range(7):
if i == 0:
nc_water_pars[i,:,:] = water_pars_array
else:
nc_water_pars[i,:,:] = water_par
water_type_array = np.ones((self.ny + 1, self.nx + 1)) * -127
for i in range(row_array_image):
for j in range(col_array_image):
if array_image[i, j] == 16.0:
water_type_array[i,j] = 1
nc_water_type = self.nc_file.createVariable(
'water_type', 'i1',('y','x'), fill_value=-127)
nc_water_type.long_name = "water type"
nc_water_type.res_orig = 1.0
#nc_water_type._FillValue = -127
nc_water_type.coordinates = "E_UTM N_UTM lon lat"
nc_water_type[:,:] = water_type_array
## resolved plant canopy
#nc_lad = self.nc_file.createVariable(
#'lad', 'f4', ('zlad','y','x'), fill_value=-9999.9)
#nc_lad.long_name = 'leaf area density'
#nc_lad[:,:,:] = -9999.9
#nc_root_area_dens_s = self.nc_file.createVariable(
#'root_area_dens_s', 'f4', ('zsoil','y','x'), fill_value=-9999.9)
#nc_root_area_dens_s.long_name = 'parameterized root area density'
#nc_root_area_dens_s.units = '1/m'
#nc_root_area_dens_s[:,:,:] = -9999.9
def finalize(self):
""" Close file """
print("Closing file...")
self.nc_file.close()
if __name__ == '__main__':
driver = StaticDriver()
driver.write_global_attributes()
driver.define_dimensions()
driver.add_variables()
driver.finalize()
以上是将地形参数写入nc文件的代码，下图显示了写入的变量以及植被的分布
最后将写好的nc文件放到项目INPUT文件夹下，重命名为MyProject_static即可。至此，静态驱动模型配置完毕。
运行结果
运行方式和自带案例的步骤一样，详情参考P2，我的运行结果有以下三个文件
下一篇讲述如何快速可视化nc文件，以及如何在arcmap中读取它。
785
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
