前言
最近需要處理一些GSMap的遙感影像，GSMaP是全球降水計(jì)劃（GPM）時(shí)代最主要的高分辨率降水產(chǎn)品，下載需要申請(qǐng)，然后通過(guò)ftp地址進(jìn)行下載。下載好以后用ENVI打開，發(fā)現(xiàn)沒有頭文件，自定義頭文件之后，影像極其的奇怪，大概是下面這個(gè)樣子，實(shí)在不符合想象，如果有同學(xué)知道怎么回事，希望您能告知我，感謝?，F(xiàn)在我打算用python先轉(zhuǎn)成tif格式。
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打開
'''
# 導(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é)果如下圖。
到此這篇關(guān)于Python讀取GSMap數(shù)據(jù)的問(wèn)題的文章就介紹到這了,更多相關(guān)Python讀取GSMap數(shù)據(jù)內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！
您可能感興趣的文章:- python讀取excel數(shù)據(jù)繪制簡(jiǎn)單曲線圖的完整步驟記錄
- Python matplotlib讀取excel數(shù)據(jù)并用for循環(huán)畫多個(gè)子圖subplot操作
- Python讀取Excel數(shù)據(jù)并生成圖表過(guò)程解析
- Python讀取JSON數(shù)據(jù)操作實(shí)例解析
- Python讀取實(shí)時(shí)數(shù)據(jù)流示例
