毫無疑問，Python是當(dāng)今最流行，最通用的編程語言之一。這有很多種強有力的原因，但在我看來，最重要的是：開源定義，語法簡單，包括電池的理念(batteries included philosophy)以及一個棒棒噠的全球社區(qū)。

Python被廣泛采用的領(lǐng)域的一個有趣的例子就是科學(xué)世界。這也解釋了像PyData 或者Scipy生態(tài)圈這類社區(qū)的存在。

另一個我看到的越來越多人感興趣并使用的更具體的區(qū)域是地理空間數(shù)據(jù)處理。這方面的一個證明是許多知名的工具使用它，例如GDAL,

ArcGIS, GRASS,QGIS等等。

因此，這篇文章的目的是為了展示在這個特定的范圍里，Python生態(tài)圈的優(yōu)勢和力量。我將通過一個復(fù)雜任務(wù)的例子來說明，這個例子是該領(lǐng)域的典型例子：衛(wèi)星圖像分類。

Python中的衛(wèi)星圖像分類

出于無數(shù)種原因，會對衛(wèi)星圖像進行分類。它用于農(nóng)業(yè)，地質(zhì)，突發(fā)事件的監(jiān)測，監(jiān)視，天氣預(yù)報，經(jīng)濟研究，社會科學(xué)等等。

出于這個原因，許多GIS和地理空間數(shù)據(jù)管理系統(tǒng)包含工具來執(zhí)行分類。但是，這種方法有一定的局限性：這個過程是手動的，你通常有一套非常小的關(guān)于分類技術(shù)及其他超參數(shù)選項。

另一種可能性是使用領(lǐng)域特定語言 (“DSL”) 進行分類，如R, IDL, MATLAB, Octave，等等。但是在這種情況下，你通常局限于一個實驗環(huán)境。

多虧了Python的腳本功能和豐富的生態(tài)系統(tǒng)，它提供了其他DSL的大多數(shù)好處，因此它非常適合于進行研究和快速原型。此外，作為一個廣泛采用的通用語言，它也有益于開發(fā)用于生產(chǎn)的，高效，可維護，可擴展，產(chǎn)業(yè)規(guī)模的分類系統(tǒng)。

因此，讓我們看看利用存在于Python生態(tài)圈中用于地理空間數(shù)據(jù)處理的工具，處理圖像分類是有多容易。在幾百行代碼中，我們將要開發(fā)的腳本：

處理一張陸地衛(wèi)星8 GeoTiff圖像(柵格數(shù)據(jù))，

從形狀文件(Shapefile)中提取訓(xùn)練和測試數(shù)據(jù)(矢量數(shù)據(jù))

使用現(xiàn)代機器學(xué)習(xí)技術(shù)訓(xùn)練和分類

評估結(jié)果

為了嘗試及kiss保持這個問題的正確分類部分的重點，我不打算深入數(shù)據(jù)預(yù)處理(校準(zhǔn)、地理變換等等)。我知道這是工作的一個基本組成部分，但它不是我現(xiàn)在要擴展的部分。

像往常一樣，大部分神奇的部分將使用作為主要依賴的工具完成：

GDAL

GDAL是一個柵格和適量地理空間數(shù)據(jù)格式的翻譯庫。對于所有支持的格式，它可以展現(xiàn)為柵格抽象數(shù)據(jù)模型或矢量抽象數(shù)據(jù)模型。

它是用C/C++實現(xiàn)的，所以具有高性能，同時，它提供了Python綁定。

安裝這個庫可能不是一個小任務(wù)，特別是對于那些不是很熟悉安裝Python依賴的過程的人來說。在任何情況下，GDAL網(wǎng)站已經(jīng)提供了詳細(xì)的指令，它們歸納在一個README文件中，位于與這篇文有關(guān)的代碼庫中。

scikit-learn

scikit-learn是一個Python開源機器學(xué)習(xí)褲。它具有多種分類、回歸和聚類算法。它被設(shè)計用于與Python數(shù)值與科學(xué)庫numpy和scipy協(xié)同工作。

它大部分使用Python編寫，一些核心算法是使用Cython(使用C/C++)編寫的，以獲得高性能。

示例數(shù)據(jù)

幸虧有了GDAL API，這篇文章中我們將開發(fā)的程序與許多不同類型的圖像格式和地理相應(yīng)的載體一起工作。但萬一你手邊沒有數(shù)據(jù)集，可以下載一些示例數(shù)據(jù)。

它包括農(nóng)業(yè)區(qū)的一張陸地衛(wèi)星8圖像，以及不同的作物樣本的某些合成矢量數(shù)據(jù)的一部分。它是那種用于精耕細(xì)作的產(chǎn)品數(shù)據(jù)(例如，作物的標(biāo)識)。

該文件是一個壓縮數(shù)據(jù)目錄，它有三個子目錄：

image 包含了一個帶有L1T陸地衛(wèi)星場景的收成的Geotiff文件(LANDSAT 8, sensor OLI, path 229, row 81, 2016-01-19 19:14:02 UTC).

只存在頻帶1到7 (Aerosol, VIS, NIR, SWIR. 30m resolution)

train 已經(jīng)帶有了用于訓(xùn)練的附帶矢量數(shù)據(jù)的shapefile。每一個文件定義了一個類別，即：所有的點，poligon等等。存在于一個shapefile中，被用于定義一個類別的樣本。在我們的樣本數(shù)據(jù)中，我們有5個列表，即A, B, C, D和E (我知道，不是太花哨。)

test類似于測試目錄，但是將使用該樣本來驗證分類結(jié)果。

由于本文將不關(guān)注分類的結(jié)果，我將不談?wù)撊魏侮P(guān)于數(shù)據(jù)質(zhì)量，要求或準(zhǔn)備的任何細(xì)節(jié)。如果你想要知道或者討論關(guān)于這個問題的任何東西，請使用評論，或者給我發(fā)郵件。

示例程序

接下來，我們將開發(fā)一個腳本來分類地理空間數(shù)據(jù)。該程序的一個更Pythonic和完備的版本可以在這個倉庫找到。那個版本包括我們在這篇文章中將不會考慮的日志，文檔字符串，評論，pep-8，一些錯誤控制和其他好的編程實踐。

在代碼庫中，你也將找到在本文中描述的腳本的一個更簡單的版本(這里)。在你下載或復(fù)制粘貼這些行到一個文件中或一個Python解析器中之前，確保安裝了下述依賴：

GDAL==2.0.1

numpy>=1.10,<1.11

scipy==0.17.0

scikit-learn==0.17

# Optionally, you can install matplotlib

準(zhǔn)備工作

現(xiàn)在，我們已經(jīng)準(zhǔn)備好寫代碼了。所以，第一件事是：導(dǎo)入我們的主要依賴，然后定義一個顏色列表：

import numpy as np

import os

from osgeo import gdal

from sklearn import metrics

from sklearn.ensemble import RandomForestClassifier

# A list of "random" colors (for a nicer output)

COLORS = ["#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941"]

從最后的導(dǎo)入行可以看到，我們將使用Random Forest技術(shù)來進行分類。多虧了Scikit-learn，我們可以輕松地實驗/比較許多不同分類技術(shù)，例如隨機梯度下降，支持向量機，最近鄰，AdaBoost等等。

顏色列表將被嵌入到GeoTiff輸出中。這將允許你輕松地在任何標(biāo)準(zhǔn)的圖像瀏覽器程序中對其可視化。

接下來，讓我們定義一些稍后將要用到的有用的函數(shù)。它們大量使用GDAL API來操縱柵格和矢量數(shù)據(jù)(代碼很好的自解釋)。

def create_mask_from_vector(vector_data_path, cols, rows, geo_transform,

projection, target_value=1):

"""Rasterize the given vector (wrapper for gdal.RasterizeLayer)."""

data_source = gdal.OpenEx(vector_data_path, gdal.OF_VECTOR)

layer = data_source.GetLayer(0)

driver = gdal.GetDriverByName('MEM') # In memory dataset

target_ds = driver.Create('', cols, rows, 1, gdal.GDT_UInt16)

target_ds.SetGeoTransform(geo_transform)

target_ds.SetProjection(projection)

gdal.RasterizeLayer(target_ds, [1], layer, burn_values=[target_value])

return target_ds

def vectors_to_raster(file_paths, rows, cols, geo_transform, projection):

"""Rasterize the vectors in the given directory in a single image."""

labeled_pixels = np.zeros((rows, cols))

for i, path in enumerate(file_paths):

label = i+1

ds = create_mask_from_vector(path, cols, rows, geo_transform,

projection, target_value=label)

band = ds.GetRasterBand(1)

labeled_pixels += band.ReadAsArray()

ds = None

return labeled_pixels

def write_geotiff(fname, data, geo_transform, projection):

"""Create a GeoTIFF file with the given data."""

driver = gdal.GetDriverByName('GTiff')

rows, cols = data.shape

dataset = driver.Create(fname, cols, rows, 1, gdal.GDT_Byte)

dataset.SetGeoTransform(geo_transform)

dataset.SetProjection(projection)

band = dataset.GetRasterBand(1)

band.WriteArray(data)

dataset = None # Close the file

現(xiàn)在，我們有了所有我們需要進行真實的分類的東西。讓我們創(chuàng)造一些變量來定義輸入和輸出：

raster_data_path = "data/image/2298119ene2016recorteTT.tif"

output_fname = "classification.tiff"

train_data_path = "data/test/"

validation_data_path = "data/train/"

在上面的幾行中，我假設(shè)我們將使用前述的樣例數(shù)據(jù)。

訓(xùn)練

現(xiàn)在，我們將使用GDAL API來讀取輸入的GeoTiff：提取地理信息，并將其轉(zhuǎn)換到一個numpy數(shù)組中：

raster_dataset = gdal.Open(raster_data_path, gdal.GA_ReadOnly)

geo_transform = raster_dataset.GetGeoTransform()

proj = raster_dataset.GetProjectionRef()

bands_data = []

for b in range(1, raster_dataset.RasterCount+1):

band = raster_dataset.GetRasterBand(b)

bands_data.append(band.ReadAsArray())

bands_data = np.dstack(bands_data)

rows, cols, n_bands = bands_data.shape

接下來，我們將處理訓(xùn)練數(shù)據(jù)：投影訓(xùn)練數(shù)據(jù)集中所有的矢量數(shù)據(jù)到一個numpy數(shù)組中。為每一個類分配一個標(biāo)簽(位于1到類的總數(shù)量之間的一個數(shù)字)。如果在這個新數(shù)組中的(i, j)位置上的值v非零，那么意味著像素(i, j)必須被用作類v的訓(xùn)練樣本。

files = [f for f in os.listdir(train_data_path) if f.endswith('.shp')]

classes = [f.split('.')[0] for f in files]

shapefiles = [os.path.join(train_data_path, f)

for f in files if f.endswith('.shp')]

labeled_pixels = vectors_to_raster(shapefiles, rows, cols, geo_transform,

proj)

is_train = np.nonzero(labeled_pixels)

training_labels = labeled_pixels[is_train]

training_samples = bands_data[is_train]

training_samples是用于訓(xùn)練的像素列表。在我們的樣例中，一個像素是帶的7維空間中的一點。

training_labels是類標(biāo)簽列表，i-th 位置表示training_samples中_i-th_像素的類。

所以現(xiàn)在，我們知道輸入圖像的哪些像素必須用于訓(xùn)練。接下來，實例化Scikit-learn中的一個RandomForestClassifier。

classifier = RandomForestClassifier(n_jobs=-1)

classifier.fit(training_samples, training_labels)

這里，有許多參數(shù)我們可以玩一玩。我建議你閱讀相關(guān)文檔，并嘗試不同的可能性。

通常情況下，這些參數(shù)的精細(xì)調(diào)諧依賴于數(shù)據(jù)，學(xué)習(xí)的特定域，記憶和處理資源，預(yù)期精度等

為了集中精力和簡單起見，我不會在這個問題上展開。正如你所看到的，我只是傳遞一個選項來使用我的電腦中的所有內(nèi)核。

分類

看！不管你信不信，這都是最難的部分。現(xiàn)在，我們有了一個訓(xùn)練模型，可以分類(預(yù)測)我們所擁有的像素數(shù)據(jù)的類別。讓我們放手去做吧。

n_samples = rows*cols

flat_pixels = bands_data.reshape((n_samples, n_bands))

result = classifier.predict(flat_pixels)

classification = result.reshape((rows, cols))

我們使用了訓(xùn)練對象來對所有的輸入圖像進行分類。我們的分類器知道如何訓(xùn)練像素，并且其predict函數(shù)期望一個像素列表，而不是一個N×M矩陣。正因為如此，我們先重塑頻帶數(shù)據(jù)，然后再分類(這使得輸出看起來像一個圖像，而不僅僅是一個多維象素的列表)。

在這一點上，如果你已經(jīng)安裝了matplotlib，你可以想像結(jié)果：

from matplotlib import pyplot as plt

f = plt.figure()

f.add_subplot(1, 2, 2)

r = bands_data[:,:,3]

g = bands_data[:,:,2]

b = bands_data[:,:,1]

rgb = np.dstack([r,g,b])

f.add_subplot(1, 2, 1)

plt.imshow(rgb/255)

f.add_subplot(1, 2, 2)

plt.imshow(classification)

但比僅僅看我們驚人的頻帶新分類更重要的是將其保存到磁盤并評估其準(zhǔn)確性。因此，讓我們做到這一點。

對于第一個任務(wù)，我們已經(jīng)創(chuàng)建了一個輔助函數(shù)：write_geotiff (你已經(jīng)忘了它，迷失在對彩色圖像的驚嘆中)。我們可以只是使用matplotlib.pyplot.imsave保存像素的數(shù)據(jù)，但之后我們會失去所有包含在GeoTiff格式中有價值的地理信息(和其他元數(shù)據(jù))。而這樣的信息對于GTS和其他衛(wèi)星數(shù)據(jù)處理系統(tǒng)是必不可少的。因此，我們將使用我們的GDAL驅(qū)動函數(shù)：

write_geotiff(output_fname, classification, geo_transform, proj)

這很簡單。現(xiàn)在，你應(yīng)該能夠打開我們剛剛創(chuàng)建的新文件，使用任何圖像瀏覽器，GIS和遙感數(shù)據(jù)處理系統(tǒng)。

評估結(jié)果

終于接近結(jié)束時，在可以驗證我們的分類的準(zhǔn)確性之前，我們需要以一種與我們處理訓(xùn)練數(shù)據(jù)相似的方式預(yù)處理我們的測試數(shù)據(jù)：

shapefiles = [os.path.join(validation_data_path, "%s.shp" % c)

for c in classes]

verification_pixels = vectors_to_raster(shapefiles, rows, cols,

geo_transform, proj)

for_verification = np.nonzero(verification_pixels)

verification_labels = verification_pixels[for_verification]

predicted_labels = classification[for_verification]

上面，我們有了要驗證像素的期望標(biāo)簽以及所計算的標(biāo)簽。因此，我們可以分析結(jié)果了。為此，我們可愛的scikit-learn提供了許多工具。所以，讓我們使用其中兩個。

print("Confussion matrix:\n%s" %

metrics.confusion_matrix(verification_labels, predicted_labels))

這應(yīng)該打印出一些像這樣的東西：

Confussion matrix:

[[ 82 0 6 0 0]

[ 0 180 0 0 0]

[ 0 0 65 0 0]

[ 0 0 2 89 0]

[ 0 0 0 0 160]]

接下來是精度和準(zhǔn)確性：

target_names = ['Class %s' % s for s in classes]

print("Classification report:\n%s" %

metrics.classification_report(verification_labels, predicted_labels,

target_names=target_names))

print("Classification accuracy: %f" %

metrics.accuracy_score(verification_labels, predicted_labels))

Should print something like this:

Classification report:

precision recall f1-score support

Class C 1.00 0.93 0.96 88

Class D 1.00 1.00 1.00 180

Class B 0.89 1.00 0.94 65

Class A 1.00 0.98 0.99 91

Class E 1.00 1.00 1.00 160

avg / total 0.99 0.99 0.99 584

Classification accuracy: 0.986301

總結(jié)

在這篇文中，我們開發(fā)了一個腳本，它處理柵格和矢量數(shù)據(jù)，使用復(fù)雜的機器學(xué)習(xí)技術(shù)進行監(jiān)督分類，可視化輸出并評估結(jié)果。

所有這一切都在一個通用的，廣泛采用的語言的100行中實現(xiàn)，并使用了高效的工具。

至少對我來說，這證明了Python生態(tài)系統(tǒng)的地理空間數(shù)據(jù)處理的好處和能力。

不幸的是，由于保密協(xié)議，我不能分享更具體的(和非常有趣的)現(xiàn)實生活中的例子。希望將來，為了擴展在這個領(lǐng)域的優(yōu)勢以及一些有用的很酷的Python技巧，我會被允許這樣做。

總結(jié)

以上是生活随笔為你收集整理的python地理数据处理 下载_python-doc/将Python用于地理空间数据处理.md at master · zhuxinyizhizun/python-doc · GitHub...的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。