快捷導(dǎo)航

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

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

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

前言

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

微信圖片_20210305235437

python庫(kù)

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

函數(shù)說(shuō)明

glob.glob:用于選擇“.dat”格式文件
np.fromfile:讀取二進(jìn)制文件到array
arcpy.NumPyArrayToRaster:將array轉(zhuǎn)化為柵格數(shù)據(jù)
arcpy.SpatialReference(4326):引用GCS_WGS_1984的空間參考代號(hào)
arcpy.DefineProjection_management:建立空間參考

完整代碼

#!/usr/bin/env python
# -*- coding:utf-8 -*-
# Name : gsmap2TIFF.py
# Author : zengsk in NanJing
# Created: 2019/8/24 23:26

'''
說(shuō)明：1.該腳本是讀取gsmap小時(shí)尺度降水?dāng)?shù)據(jù)，輸出為tiff
 2.運(yùn)行環(huán)境需要安裝python2 需要arcpy模塊
 3.使用arcgis自帶的python環(huán)境(有arcpy模塊)
 4.運(yùn)行結(jié)果可以直接用arcgis打開(kāi)
'''

# 導(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

 # 輸出為T(mén)IFF(注意：要用到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é)果如下圖。

到此這篇關(guān)于Python讀取GSMap數(shù)據(jù)的問(wèn)題的文章就介紹到這了,更多相關(guān)Python讀取GSMap數(shù)據(jù)內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！

您可能感興趣的文章: