之前用idl寫高分預處理的時候，就有想過可不可以用python+gdal寫，可是一直卡在了第一步的正射校正，gdal.warp()函式始終找不到放dem的位置，最近終於找到了。我嘗試了一景1.3g的gf1/wfv，採用envi/idl的指令碼執行每次都需要500s以上，而python3+osgeo則穩定在驚人的15s以內！就速度而言，python3+osgeo遠遠快於envi介面。
以下是今天寫的簡單的**，包括解壓函式，正射校正函式和融合函式（gdal的融合方法只有預設的加權brovey變換）。執行了一景gf6/pms，並和idl版白鴿的結果作了簡單對比。
import sys, os, tarfile, tempfile as tmp
from osgeo import gdal, osr
from gdalconst import
*def
unpackage
(fn)
: dirname, basename = os.path.split(fn)
odir = os.path.join(dirname,
'snowydove_'
+ basename[0:
-7])
with tarfile.
open
(fn)
asfile
:file
.extractall(path = odir)
return odir
defortho
(ifn, demfn, res, ofn)
: ds = gdal.open(ifn, ga_readonly)
isnorth =
1if os.path.basename(ifn)
.split(
'_')[3
][0]
=='n'
else
0 zone =
str(
int(
float
(os.path.basename(ifn)
.split(
'_')[2
][1:
])/6
)+31)
zone =
int(
'326'
+ zone)
if isnorth else
int(
'327'
+ zone)
dstsrs = osr.spatialreference(
) dstsrs.importfromepsg(zone)
tds = gdal.warp(ofn, ds,
format
='gtiff'
, xres = res, yres = res, dstsrs = dstsrs, rpc =
true
, transformeroptions = demfn)
ds = tds =
none
deffindimage
(dir):
fn =
with os.scandir(
dir)
as it:
for entry in it:
if entry.name.endswith(
'.tiff'):
dir, entry.name))if
len(fn)==1
:return fn[0]
elif
len(fn)==2
:if fn[0]
.find(
'pan')!=
-1:return fn[1]
, fn[0]
else
:return fn[0]
, fn[1]
elif
len(fn)==3
:return fn
elif
len(fn)==6
: mss = pan =
none
for element in fn:
if element.find(
'mux.')!=
-1: mss = element
continue
if element.find(
'pan.')!=
-1: pan = element
continue
if mss !=
none
and pan !=
none
:break
return mss, pan
defpansharpen
(mss, pan, ofn)
: ds = gdal.open(mss, ga_readonly)
nb = ds.rastercount
ds =
none
vrt =
"""%s
1\n"""
% pan
for i in
range
(nb)
: vrt +=
""" %s
%s\n"""
%(i+
1, mss, i+1)
vrt +=
"""
\n"""
pansharpends = gdal.open(vrt)
newds = gdal.getdriverbyname(
'gtiff'
).createcopy(ofn, pansharpends)
pansharpends = newds =
none
defsnowydove
(ar**)
: imgdir = unpackage(sys.ar**[1]
) mss, pan = findimage(imgdir)
odir = os.path.join(os.path.dirname(sys.ar**[1]
),'snowydove_tempfiles'
) os.mkdir(odir)
ortho_mss_fn = tmp.mkstemp(
'.tiff'
,'snowydove_'
, odir,
false)[
1]ortho(mss, sys.ar**[2]
,8, ortho_mss_fn)
print
('mss done'
) ortho_pan_fn = tmp.mkstemp(
'.tiff'
,'snowydove_'
, odir,
false)[
1]ortho(pan, sys.ar**[2]
,2, ortho_pan_fn)
print
('pan done'
) sharpen_fn = tmp.mkstemp(
'.tiff'
,'snowydove_'
, odir,
false)[
1]pansharpen(ortho_mss_fn, ortho_pan_fn, sharpen_fn)
return
0def
main()
:return snowydove(sys.ar**)
if __name__ ==
'__main__'
: sys.exit(snowydove(sys.ar**)
)
後來分析才發現，其實產生這種錯位差異的原因可能在於正射校正環節：gdal的正射校正保持了原汁原味的鋸齒狀，而envi的正射校正則使用了一定的平滑，所以才可能導致全色和多光譜的錯位。
Python處理高光譜資料 3 降維
用到的庫 matplotlib scipy spectral 主要內容 主成分分析 pca 與線性判別 lda 歡迎有興趣的朋友交流指點。最後，廢話不多說直接上 import matplotlib.pyplot as plt from scipy.io import loadmat import s...
python 異常處理3
def set age age if age 0 or age 200 print 值錯誤 raise valueerror 值錯誤 else print 給張三的年齡設定為 age try set age 18 except exception as e print x e 什麼時候應該向外界丟擲...
python學習高光時刻記錄3
1.title 將每個單詞的首字母都改為大寫 upper 將字串改為全部大寫 lower 將字串改為全部小寫 title name love you print name.title 輸出 love you upper name love you print name.upper 輸出 love y...