爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 1025|回复: 2

[求助] 三维的降水数据如何进行区域的选取

[复制链接]

新浪微博达人勋

发表于 2023-12-26 14:39:53 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
我想要对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. # 读取TRMM的HDF4文件
  7. file_path = '/home/19980101.hdf'

  8. # 打开HDF4文件并获取数据集
  9. hdf_file = SD(file_path, SDC.READ)

  10. # 获取经纬度数据
  11. longitude_data = np.array(hdf_file.select('Longitude').get())
  12. latitude_data = np.array(hdf_file.select('Latitude').get())

  13. # 假设你已经准备好了包含省份边界信息的 GeoJSON 文件
  14. geojson_file_path = '/home/四川省.json'

  15. # 读取省份边界数据
  16. with open(geojson_file_path) as f:
  17.     province_data = json.load(f)

  18. # 提取省份边界数据的polygon
  19. province_polygon = shape(province_data['features'][0]['geometry'])

  20. # 获取降水数据
  21. rain_data = hdf_file.select('rain').get()

  22. # 获取 rain 数据的形状
  23. rain_shape = rain_data.shape

  24. # 创建一个列表来存储每层高度的保留数据
  25. filtered_data_by_height = []

  26. # 遍历每个高度层
  27. for h in range(rain_shape[2]):
  28.     points_in_province = []

  29.     # 遍历经度和纬度数据
  30.     for i in range(rain_shape[0]):
  31.         for j in range(rain_shape[1]):
  32.             point = Point(longitude_data[i][j], latitude_data[i][j])
  33.             if province_polygon.contains(point):
  34.                 points_in_province.append(rain_data[i][j][h])

  35.     # 将符合条件的降水数据转换为 Numpy 数组并存储
  36.     filtered_data_by_height.append(np.array(points_in_province))
  37.    
  38. # 使用 np.dstack() 合并每个高度层的保留数据
  39. combined_data = np.dstack(filtered_data_by_height)

  40. # 确定输出文件路径和文件名
  41. output_directory = '/home/19980101裁剪区域是四川.hdf'

  42. # 获取目录路径
  43. output_dir = os.path.dirname(output_directory)

  44. # 如果目录不存在,则创建
  45. if not os.path.exists(output_dir):
  46.     os.makedirs(output_dir)

  47. # 创建新的HDF4文件
  48. output_hdf = SD(output_directory, SDC.WRITE | SDC.CREATE)

  49. # 将数据写入新文件
  50. for idx, layer_data in enumerate(filtered_data_by_height):
  51.     output_hdf.create(f'layer_{idx+1}', SDC.FLOAT64, layer_data.shape)
  52.     output_hdf.select(f'layer_{idx+1}').set(layer_data)

  53. # 关闭文件
  54. output_hdf.end()

  55. print('successful')
复制代码


密码修改失败请联系微信:mofangbao

新浪微博达人勋

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

不过这个教程比较老,现在用pyproj、geopadas做mask更方便
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

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

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

我按您发的贴子和思路试一下!感谢!~
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表