三维的降水数据如何进行区域的选取

小k

我想要对TRMM 2A25的二级数据（它是一个三维的降水数据，包括经纬度和高度信息）的范围做一个裁剪，只留下四川地区的数据，目前已经准备好了json数据什么的，但是由于这个数据是三维的，它的高度是一层一层的，所以不太会弄。/(ㄒoㄒ)/~~

写了一个代码，但是运行的不太好，得出来的裁剪好的数据变成了二维的😵

不知道哪里写的不对，想请各位大佬们帮忙看一下，感谢感谢！

1. from pyhdf.SD import SD, SDC
   
2. from shapely.geometry import shape, Point
   
3. import json
   
4. import numpy as np
   
5. import os
   
6. 
   
7. # 读取TRMM的HDF4文件
   
8. file_path = '/home/19980101.hdf'
   
9. 
   
10. # 打开HDF4文件并获取数据集
    
11. hdf_file = SD(file_path, SDC.READ)
    
12. 
    
13. # 获取经纬度数据
    
14. longitude_data = np.array(hdf_file.select('Longitude').get())
    
15. latitude_data = np.array(hdf_file.select('Latitude').get())
    
16. 
    
17. # 假设你已经准备好了包含省份边界信息的 GeoJSON 文件
    
18. geojson_file_path = '/home/四川省.json'
    
19. 
    
20. # 读取省份边界数据
    
21. with open(geojson_file_path) as f:
    
22.     province_data = json.load(f)
    
23. 
    
24. # 提取省份边界数据的polygon
    
25. province_polygon = shape(province_data['features'][0]['geometry'])
    
26. 
    
27. # 获取降水数据
    
28. rain_data = hdf_file.select('rain').get()
    
29. 
    
30. # 获取 rain 数据的形状
    
31. rain_shape = rain_data.shape
    
32. 
    
33. # 创建一个列表来存储每层高度的保留数据
    
34. filtered_data_by_height = []
    
35. 
    
36. # 遍历每个高度层
    
37. for h in range(rain_shape[2]):
    
38.     points_in_province = []
    
39. 
    
40.     # 遍历经度和纬度数据
    
41.     for i in range(rain_shape[0]):
    
42.         for j in range(rain_shape[1]):
    
43.             point = Point(longitude_data[i][j], latitude_data[i][j])
    
44.             if province_polygon.contains(point):
    
45.                 points_in_province.append(rain_data[i][j][h])
    
46. 
    
47.     # 将符合条件的降水数据转换为 Numpy 数组并存储
    
48.     filtered_data_by_height.append(np.array(points_in_province))
    
49.    
    
50. # 使用 np.dstack() 合并每个高度层的保留数据
    
51. combined_data = np.dstack(filtered_data_by_height)
    
52. 
    
53. # 确定输出文件路径和文件名
    
54. output_directory = '/home/19980101裁剪区域是四川.hdf'
    
55. 
    
56. # 获取目录路径
    
57. output_dir = os.path.dirname(output_directory)
    
58. 
    
59. # 如果目录不存在，则创建
    
60. if not os.path.exists(output_dir):
    
61.     os.makedirs(output_dir)
    
62. 
    
63. # 创建新的HDF4文件
    
64. output_hdf = SD(output_directory, SDC.WRITE | SDC.CREATE)
    
65. 
    
66. # 将数据写入新文件
    
67. for idx, layer_data in enumerate(filtered_data_by_height):
    
68.     output_hdf.create(f'layer_{idx+1}', SDC.FLOAT64, layer_data.shape)
    
69.     output_hdf.select(f'layer_{idx+1}').set(layer_data)
    
70. 
    
71. # 关闭文件
    
72. output_hdf.end()
    
73. 
    
74. print('successful')
    

复制代码

墨家大宝

每一层都重复判断该格点在不在研究区内也太费时了吧，最后得到的而二维数组每一行是一个高度层，水平位置信息全部丢失。建议构建一个和你数据匹配的研究区域mask，参考http://bbs.06climate.com/forum.php?mod=viewthread&tid=96376

不过这个教程比较老，现在用pyproj、geopadas做mask更方便

小k

墨家大宝 发表于 2023-12-26 17:22
每一层都重复判断该格点在不在研究区内也太费时了吧，最后得到的而二维数组每一行是一个高度层，水平位置信 ...

是的您说的对！它的代码跑得很慢，最后得到一个二位数组，我之前一直没理解/(ㄒoㄒ)/~~

我按您发的贴子和思路试一下！感谢！~