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

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

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

ArcGIS, GRASS,QGIS等等。

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

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

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

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

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

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

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

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

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

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

評(píng)估結(jié)果

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

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

GDAL

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

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

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

scikit-learn

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

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

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

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

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

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

image 包含了一個(gè)帶有L1T陸地衛(wèi)星場(chǎng)景的收成的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。每一個(gè)文件定義了一個(gè)類(lèi)別，即：所有的點(diǎn)，poligon等等。存在于一個(gè)shapefile中，被用于定義一個(gè)類(lèi)別的樣本。在我們的樣本數(shù)據(jù)中，我們有5個(gè)列表，即A, B, C, D和E (我知道，不是太花哨。)

test類(lèi)似于測(cè)試目錄，但是將使用該樣本來(lái)驗(yàn)證分類(lèi)結(jié)果。

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

示例程序

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

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

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)備好寫(xiě)代碼了。所以，第一件事是：導(dǎo)入我們的主要依賴(lài)，然后定義一個(gè)顏色列表：

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ù)來(lái)進(jìn)行分類(lèi)。多虧了Scikit-learn，我們可以輕松地實(shí)驗(yàn)/比較許多不同分類(lèi)技術(shù)，例如隨機(jī)梯度下降，支持向量機(jī)，最近鄰，AdaBoost等等。

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

接下來(lái)，讓我們定義一些稍后將要用到的有用的函數(shù)。它們大量使用GDAL API來(lái)操縱柵格和矢量數(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)在，我們有了所有我們需要進(jìn)行真實(shí)的分類(lèi)的東西。讓我們創(chuàng)造一些變量來(lái)定義輸入和輸出：

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來(lái)讀取輸入的GeoTiff：提取地理信息，并將其轉(zhuǎn)換到一個(gè)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

接下來(lái)，我們將處理訓(xùn)練數(shù)據(jù)：投影訓(xùn)練數(shù)據(jù)集中所有的矢量數(shù)據(jù)到一個(gè)numpy數(shù)組中。為每一個(gè)類(lèi)分配一個(gè)標(biāo)簽(位于1到類(lèi)的總數(shù)量之間的一個(gè)數(shù)字)。如果在這個(gè)新數(shù)組中的(i, j)位置上的值v非零，那么意味著像素(i, j)必須被用作類(lèi)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)練的像素列表。在我們的樣例中，一個(gè)像素是帶的7維空間中的一點(diǎn)。

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

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

classifier = RandomForestClassifier(n_jobs=-1)

classifier.fit(training_samples, training_labels)

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

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

為了集中精力和簡(jiǎn)單起見(jiàn)，我不會(huì)在這個(gè)問(wèn)題上展開(kāi)。正如你所看到的，我只是傳遞一個(gè)選項(xiàng)來(lái)使用我的電腦中的所有內(nèi)核。

分類(lèi)

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

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)練對(duì)象來(lái)對(duì)所有的輸入圖像進(jìn)行分類(lèi)。我們的分類(lèi)器知道如何訓(xùn)練像素，并且其predict函數(shù)期望一個(gè)像素列表，而不是一個(gè)N×M矩陣。正因?yàn)槿绱?#xff0c;我們先重塑頻帶數(shù)據(jù)，然后再分類(lèi)(這使得輸出看起來(lái)像一個(gè)圖像，而不僅僅是一個(gè)多維象素的列表)。

在這一點(diǎn)上，如果你已經(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)

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

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

write_geotiff(output_fname, classification, geo_transform, proj)

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

評(píng)估結(jié)果

終于接近結(jié)束時(shí)，在可以驗(yàn)證我們的分類(lèi)的準(zhǔn)確性之前，我們需要以一種與我們處理訓(xùn)練數(shù)據(jù)相似的方式預(yù)處理我們的測(cè)試數(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]

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

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

接下來(lái)是精度和準(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é)

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

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

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

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

總結(jié)

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

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