Python读取GSMap数据的问题
前言
最近需要处理一些GSMap的遥感影像，GSMaP是全球降水计划（GPM）时代最主要的高分辨率降水产品，下载需要申请，然后通过ftp地址进行下载。下载好以后用ENVI打开，发现没有头文件，自定义头文件之后，影像极其的奇怪，大概是下面这个样子，实在不符合想象，如果有同学知道怎么回事，希望您能告知我，感谢。现在我打算用python先转成tif格式。
python库
glob:是一个文件操作相关模块，支持通配符操作，可用于查找符合一定条件的文件，这里用来批量查找".dat"文件 numpy:读取“.dat”文件为array，固定尺寸为[1200,3600]，填补异常值和缺失值为-999.0 arcpy:将numpy的array输出为“.tiff”格式，添加空间参考GCS_WGS_1984，最终保存tif文件
函数说明
glob.glob:用于选择“.dat”格式文件 np.fromfile:读取二进制文件到array arcpy.NumPyArrayToRaster:将array转化为栅格数据 arcpy.SpatialReference(4326):引用GCS_WGS_1984的空间参考代号 arcpy.DefineProjection_management:建立空间参考
完整代码
#!/usr/bin/envpython
#-*-coding:utf-8-*-
#Name:gsmap2TIFF.py
#Author:zengskinNanJing
#Created:2019/8/2423:26
'''
说明：1.该脚本是读取gsmap小时尺度降水数据，输出为tiff
2.运行环境需要安装python2需要arcpy模块
3.使用arcgis自带的python环境(有arcpy模块)
4.运行结果可以直接用arcgis打开
'''
#导入module
importos
importglob
importnumpyasnp
importarcpy
importwarnings
warnings.simplefilter("ignore")#忽略警告
#原始降水数据文件夹,可根据自己本地情况修改
sPath=r'./data'
oDir=r"./tiff/"
forfileNameinglob.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模块)
ifnotos.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++++++DataProcessingSuccessfullyCompleted!++++++")
最终结果
最终得到GSMap的栅格数据，示例结果如下图。
到此这篇关于Python读取GSMap数据的问题的文章就介绍到这了,更多相关Python读取GSMap数据内容请搜索毛票票以前的文章或继续浏览下面的相关文章希望大家以后多多支持毛票票！
声明：本文内容来源于网络，版权归原作者所有，内容由互联网用户自发贡献自行上传，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任。如果您发现有涉嫌版权的内容，欢迎发送邮件至：czq8825#qq.com（发邮件时，请将#更换为@）进行举报，并提供相关证据，一经查实，本站将立刻删除涉嫌侵权内容。