登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 晋陵小生 于 2022-2-23 22:55 编辑
本文已发布于微信公众号【气象备忘录】,所有代码和数据均同步至和鲸社区
传送门:公众号文章;在线代码
旧文新发,1. 修复了数据偏移的问题;2. 操作更加简单,fork代码后只需要修改经纬度就能在线制作WPS可用的地形数据。以下是正文部分。
利用ASTER-30米地形数据,制作WPS中可以使用的高分辨率地形高度数据,可以为高分辨率数值模拟(如大涡)提供较为准确的地形高度信息。 某度网盘下载速度实在感人,我直接把数据搬到和鲸了。大家可以自行设定经纬度,在线完成数据的制作,最后打包下载到本地使用。上传的数据集中包括中国区域的数据,只需要设定好经纬度范围即可,其他内容无需修改。 处理思路主要参考: http://bbs.06climate.com/forum.php?mod=viewthread&tid=50088&highlight=aster
tips1: 上面链接存在一个问题,两个数据块之间存在一条重复记录。即每拼接一块数据就会产生一条记录的偏移,因此需要在格式转换的时候剔除重复记录(感谢宝哥@内蒙气象局)。
tips2: 选取范围不能过大,静态数据命名上限为99999,1个数据块3600,即范围不能超过27个经纬度 。
以下是全部代码: - import os
- import shutil
- import zipfile
- import glob
- # 需要解压的文件夹
- ZIP_ROOT = "/home/mw/input"
- # 解压到的根目录
- EXTRACT_BASE_PATH = r"/home/mw/project/Extracted"
- aster_dir = "/home/mw/project/topo_ASTER"
- if not os.path.exists(aster_dir):
- os.makedirs(aster_dir)
- if not os.path.exists(EXTRACT_BASE_PATH):
- os.makedirs(EXTRACT_BASE_PATH)
- lat1, lat2, lon1, lon2 = 40, 42, 115, 117
- ziplist = [ f"ASTGTMV003_N{lat}E{lon}.zip" for lat in range(lat1, lat2+1) for lon in range(lon1, lon2+1) ]
- for root, dirs, files in os.walk(ZIP_ROOT):
- for file in files:
- if file in ziplist:
- fn = os.path.join(root,file)
- print(f"正在解压->{fn}")
- zip = zipfile.ZipFile(fn,'r')
- zip.extractall(path=EXTRACT_BASE_PATH)
- zip.close()
- os.environ['PROJ_LIB'] = "/opt/conda/share/proj" # 不加这句会报proj的错
- files = glob.glob(f"{EXTRACT_BASE_PATH}/*dem.tif")
- files.sort()
- for infile in files:
- print(infile)
- outfile = infile[:-3] + "bil"
- cmd = f"/opt/conda/bin/gdal_translate -of ENVI -co INTERLEAVE=BSQ -srcwin 1 0 3600 3600 {infile} {outfile}"
- # print(cmd)
- os.system(cmd)
- files = glob.glob(f"{EXTRACT_BASE_PATH}/*bil")
- files.sort()
- for ifile in files:
- lat = int(ifile[-14:-12])
- lon = int(ifile[-11:-8])
- x1 = (lon - lon1)*3600 + 1
- x2 = (lon - lon1 + 1)*3600
- y1 = (lat - lat1)*3600 + 1
- y2 = (lat - lat1 + 1)*3600
- to_ = os.path.join(aster_dir, f"{x1:05d}-{x2:05d}.{y1:05d}-{y2:05d}")
- shutil.move(ifile,to_ )
- index = f'''\
- type = continuous
- signed = yes
- projection = regular_ll
- dx = 0.00027777778
- dy = 0.00027777778
- known_x = 1.0
- known_y = 1.0
- known_lat = {lat1}.000
- known_lon = {lon1}.000
- wordsize = 2
- endian = little
- tile_x = 3600
- tile_y = 3600
- tile_z = 1
- row_order = top_bottom
- missing_value = 32768
- units = "meters MSL"
- description = "ASTER 1-sec Topography Height"
- '''
- with open(os.path.join(aster_dir, 'index'),'w') as f:
- f.write(index)
- # cd /home/mw/project && tar -cvf aster.tar topo_ASTER
复制代码欢迎fork&测试
|