快捷導(dǎo)航

Python讀取GSMap數(shù)據(jù)的問題

更新時間：2021年03月10日 10:17:41 作者：凡人不修仙

這篇文章主要介紹了Python讀取GSMap數(shù)據(jù)的問題,本文給大家介紹的非常詳細(xì)，對大家的學(xué)習(xí)或工作具有一定的參考借鑒價值，需要的朋友可以參考下

前言

最近需要處理一些GSMap的遙感影像，GSMaP是全球降水計劃（GPM）時代最主要的高分辨率降水產(chǎn)品，下載需要申請，然后通過ftp地址進(jìn)行下載。下載好以后用ENVI打開，發(fā)現(xiàn)沒有頭文件，自定義頭文件之后，影像極其的奇怪，大概是下面這個樣子，實(shí)在不符合想象，如果有同學(xué)知道怎么回事，希望您能告知我，感謝?，F(xiàn)在我打算用python先轉(zhuǎn)成tif格式。

微信圖片_20210305235437

python庫

glob:是一個文件操作相關(guān)模塊，支持通配符操作，可用于查找符合一定條件的文件，這里用來批量查找".dat"文件
numpy:讀取“.dat”文件為array，固定尺寸為[1200, 3600]，填補(bǔ)異常值和缺失值為-999.0
arcpy:將numpy的array輸出為“.tiff”格式，添加空間參考GCS_WGS_1984，最終保存tif文件

函數(shù)說明

glob.glob:用于選擇“.dat”格式文件
np.fromfile:讀取二進(jìn)制文件到array
arcpy.NumPyArrayToRaster:將array轉(zhuǎn)化為柵格數(shù)據(jù)
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小時尺度降水?dāng)?shù)據(jù)，輸出為tiff
 2.運(yùn)行環(huán)境需要安裝python2 需要arcpy模塊
 3.使用arcgis自帶的python環(huán)境(有arcpy模塊)
 4.運(yùn)行結(jié)果可以直接用arcgis打開
'''

# 導(dǎo)入module
import os
import glob
import numpy as np
import arcpy
import warnings

warnings.simplefilter("ignore") # 忽略警告

# 原始降水?dāng)?shù)據(jù)文件夾,可根據(jù)自己本地情況修改
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' # 輸出文件名(可根據(jù)實(shí)際情況改)
 # arcpy.NumPyArrayToRaster()不清楚輸入?yún)?shù)可以查看arcpy的官方文檔
 # 矩陣轉(zhuǎn)為柵格
 raster = arcpy.NumPyArrayToRaster(ds, arcpy.Point(0, -60.0),
          x_cell_size=0.1, y_cell_size=0.1, value_to_nodata=-999.00)
 # 添加地理坐標(biāo)系 GCS_WGS_1984
 spatialRef = arcpy.SpatialReference(4326)
 arcpy.DefineProjection_management(raster, spatialRef)
 raster.save(TiffName)

print("\n++++++ Data Processing Successfully Completed ! ++++++")

最終結(jié)果

最終得到GSMap的柵格數(shù)據(jù)，示例結(jié)果如下圖。