Python读取GSMap数据的问题-6ren

Python读取GSMap数据的问题

转载作者：qq735679552 更新时间：2022-09-29 22:32:09

CFSDN坚持开源创造价值，我们致力于搭建一个资源共享平台，让每一个IT人在这里找到属于你的精彩世界.

这篇CFSDN的博客文章Python读取GSMap数据的问题由作者收集整理，如果你对这篇文章有兴趣，记得点赞哟.

前言

最近需要处理一些GSMap的遥感影像，GSMaP是全球降水计划（GPM）时代最主要的高分辨率降水产品，下载需要申请，然后通过ftp地址进行下载。下载好以后用ENVI打开，发现没有头文件，自定义头文件之后，影像极其的奇怪，大概是下面这个样子，实在不符合想象，如果有同学知道怎么回事，希望您能告知我，感谢。现在我打算用python先转成tif格式.

Python读取GSMap数据的问题

python库

glob:是一个文件操作相关模块，支持通配符操作，可用于查找符合一定条件的文件，这里用来批量查找".dat"文件numpy:读取“.dat”文件为array，固定尺寸为[1200, 3600]，填补异常值和缺失值为-999.0arcpy:将numpy的array输出为“.tiff”格式，添加空间参考GCS_WGS_1984，最终保存tif文件

函数说明

glob.glob:用于选择“.dat”格式文件np.fromfile:读取二进制文件到arrayarcpy.NumPyArrayToRaster:将array转化为栅格数据arcpy.SpatialReference(4326):引用GCS_WGS_1984的空间参考代号arcpy.DefineProjection_management:建立空间参考

完整代码。

#!/usr/bin/env python# -*- coding:utf-8 -*-# Name : gsmap2TIFF.py# Author : zengsk in NanJing# Created: 2019/8/24 23:26'''说明：1.该脚本是读取gsmap小时尺度降水数据，输出为tiff 2.运行环境需要安装python2 需要arcpy模块 3.使用arcgis自带的python环境(有arcpy模块) 4.运行结果可以直接用arcgis打开'''# 导入moduleimport osimport globimport numpy as npimport arcpyimport warningswarnings.simplefilter("ignore") # 忽略警告# 原始降水数据文件夹,可根据自己本地情况修改sPath = r'./data'oDir = r"./tiff/"for fileName in glob.glob(sPath + '\*.dat'): print("Processing... {0}".format(fileName)) ds = np.fromfile(fileName, dtype=np.float32) ds = np.resize(ds, [1200, 3600]) ds[ds < 0] = -999.00 ds[np.isnan(ds)] = -999.00 # NODATA_value # 输出为TIFF(注意：要用到arcpy模块) if not os.path.exists(oDir):  os.makedirs(oDir) TiffName = oDir + os.sep + os.path.basename(fileName)[0:-4] + '.tif' # 输出文件名(可根据实际情况改) # arcpy.NumPyArrayToRaster()不清楚输入参数可以查看arcpy的官方文档 # 矩阵转为栅格 raster = arcpy.NumPyArrayToRaster(ds, arcpy.Point(0, -60.0),          x_cell_size=0.1, y_cell_size=0.1, value_to_nodata=-999.00) # 添加地理坐标系 GCS_WGS_1984 spatialRef = arcpy.SpatialReference(4326) arcpy.DefineProjection_management(raster, spatialRef) raster.save(TiffName)print("\n++++++ Data Processing Successfully Completed ! ++++++")

最终结果。

最终得到GSMap的栅格数据，示例结果如下图.

Python读取GSMap数据的问题