毫無疑問，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技巧，我會被允許這樣做。