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

Python被廣泛采用的領域的一個有趣的例子就是科學世界。這也解釋了像PyData 或者Scipy生態圈這類社區的存在。

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

ArcGIS, GRASS,QGIS等等。

因此，這篇文章的目的是為了展示在這個特定的范圍里，Python生態圈的優勢和力量。我將通過一個復雜任務的例子來說明，這個例子是該領域的典型例子：衛星圖像分類。

Python中的衛星圖像分類

出于無數種原因，會對衛星圖像進行分類。它用于農業，地質，突發事件的監測，監視，天氣預報，經濟研究，社會科學等等。

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

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

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

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

處理一張陸地衛星8 GeoTiff圖像(柵格數據)，

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

使用現代機器學習技術訓練和分類

評估結果

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

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

GDAL

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

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

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

scikit-learn

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

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

示例數據

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

它包括農業區的一張陸地衛星8圖像，以及不同的作物樣本的某些合成矢量數據的一部分。它是那種用于精耕細作的產品數據(例如，作物的標識)。

該文件是一個壓縮數據目錄，它有三個子目錄：

image 包含了一個帶有L1T陸地衛星場景的收成的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 已經帶有了用于訓練的附帶矢量數據的shapefile。每一個文件定義了一個類別，即：所有的點，poligon等等。存在于一個shapefile中，被用于定義一個類別的樣本。在我們的樣本數據中，我們有5個列表，即A, B, C, D和E (我知道，不是太花哨。)

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

由于本文將不關注分類的結果，我將不談論任何關于數據質量，要求或準備的任何細節。如果你想要知道或者討論關于這個問題的任何東西，請使用評論，或者給我發郵件。

示例程序

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

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

GDAL==2.0.1

numpy>=1.10,<1.11

scipy==0.17.0

scikit-learn==0.17

# Optionally, you can install matplotlib

準備工作

現在，我們已經準備好寫代碼了。所以，第一件事是：導入我們的主要依賴，然后定義一個顏色列表：

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"]

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

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

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

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

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

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

output_fname = "classification.tiff"

train_data_path = "data/test/"

validation_data_path = "data/train/"

在上面的幾行中，我假設我們將使用前述的樣例數據。

訓練

現在，我們將使用GDAL API來讀取輸入的GeoTiff：提取地理信息，并將其轉換到一個numpy數組中：

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

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

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是用于訓練的像素列表。在我們的樣例中，一個像素是帶的7維空間中的一點。

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

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

classifier = RandomForestClassifier(n_jobs=-1)

classifier.fit(training_samples, training_labels)

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

通常情況下，這些參數的精細調諧依賴于數據，學習的特定域，記憶和處理資源，預期精度等

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

分類

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

n_samples = rows*cols

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

result = classifier.predict(flat_pixels)

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

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

在這一點上，如果你已經安裝了matplotlib，你可以想像結果：

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)

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

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

write_geotiff(output_fname, classification, geo_transform, proj)

這很簡單。現在，你應該能夠打開我們剛剛創建的新文件，使用任何圖像瀏覽器，GIS和遙感數據處理系統。

評估結果

終于接近結束時，在可以驗證我們的分類的準確性之前，我們需要以一種與我們處理訓練數據相似的方式預處理我們的測試數據：

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]

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

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

metrics.confusion_matrix(verification_labels, predicted_labels))

這應該打印出一些像這樣的東西：

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]]

接下來是精度和準確性：

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

總結

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

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

至少對我來說，這證明了Python生態系統的地理空間數據處理的好處和能力。

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

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。