- 积分
- 30
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2023-6-2
- 最后登录
- 1970-1-1
![未绑定新浪微博用户 新浪微博达人勋](source/plugin/sina_login/img/gray.png)
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
我想要对TRMM 2A25的二级数据(它是一个三维的降水数据,包括经纬度和高度信息)的范围做一个裁剪,只留下四川地区的数据,目前已经准备好了json数据什么的,但是由于这个数据是三维的,它的高度是一层一层的,所以不太会弄。/(ㄒoㄒ)/~~
写了一个代码,但是运行的不太好,得出来的裁剪好的数据变成了二维的😵
不知道哪里写的不对,想请各位大佬们帮忙看一下,感谢感谢!
- from pyhdf.SD import SD, SDC
- from shapely.geometry import shape, Point
- import json
- import numpy as np
- import os
- # 读取TRMM的HDF4文件
- file_path = '/home/19980101.hdf'
- # 打开HDF4文件并获取数据集
- hdf_file = SD(file_path, SDC.READ)
- # 获取经纬度数据
- longitude_data = np.array(hdf_file.select('Longitude').get())
- latitude_data = np.array(hdf_file.select('Latitude').get())
- # 假设你已经准备好了包含省份边界信息的 GeoJSON 文件
- geojson_file_path = '/home/四川省.json'
- # 读取省份边界数据
- with open(geojson_file_path) as f:
- province_data = json.load(f)
- # 提取省份边界数据的polygon
- province_polygon = shape(province_data['features'][0]['geometry'])
- # 获取降水数据
- rain_data = hdf_file.select('rain').get()
- # 获取 rain 数据的形状
- rain_shape = rain_data.shape
- # 创建一个列表来存储每层高度的保留数据
- filtered_data_by_height = []
- # 遍历每个高度层
- for h in range(rain_shape[2]):
- points_in_province = []
- # 遍历经度和纬度数据
- for i in range(rain_shape[0]):
- for j in range(rain_shape[1]):
- point = Point(longitude_data[i][j], latitude_data[i][j])
- if province_polygon.contains(point):
- points_in_province.append(rain_data[i][j][h])
- # 将符合条件的降水数据转换为 Numpy 数组并存储
- filtered_data_by_height.append(np.array(points_in_province))
-
- # 使用 np.dstack() 合并每个高度层的保留数据
- combined_data = np.dstack(filtered_data_by_height)
- # 确定输出文件路径和文件名
- output_directory = '/home/19980101裁剪区域是四川.hdf'
- # 获取目录路径
- output_dir = os.path.dirname(output_directory)
- # 如果目录不存在,则创建
- if not os.path.exists(output_dir):
- os.makedirs(output_dir)
- # 创建新的HDF4文件
- output_hdf = SD(output_directory, SDC.WRITE | SDC.CREATE)
- # 将数据写入新文件
- for idx, layer_data in enumerate(filtered_data_by_height):
- output_hdf.create(f'layer_{idx+1}', SDC.FLOAT64, layer_data.shape)
- output_hdf.select(f'layer_{idx+1}').set(layer_data)
- # 关闭文件
- output_hdf.end()
- print('successful')
复制代码
|
|