什麼是柵格數據？

最簡形式的柵格由按行和列（或格網）組織的像元（或畫素）矩陣組成，其中的每個像元都包含一個資訊值（例如溫度）。柵格可以是數位航空像片、衛星影像、數位圖片或甚至掃描的地圖。

爲何將數據儲存爲柵格？

在GIS的應用中最常見的是向量數據和柵格數據，相比於向量數據，柵格數據的儲存格式簡單，處理簡單，所以經常使用。但是也存在着數據冗餘大的缺點。

有時只能將數據儲存爲柵格；
例如，影像僅以柵格形式提供。然而，許多其他要素（例如點要素）和測量值（例如降雨量）既可以儲存爲柵格數據型別也可以儲存爲要素（向量）數據型別。

柵格數據的一般特徵

柵格數據集中，每個像元都有一個值，用來表達是所描繪的現象，如類別、高度、量級或光譜等等。

在柵格數據集中，每個像元（也稱爲畫素）都有一個值。此像元值表示的是柵格數據集所描繪的現象，如類別、量級、高度或光譜值等。而其中的類別則可以是草地、森林或道路等土地利用類。量級可以表示重力、噪聲污染或降雨百分比。高度（距離）則可表示平均海平面以上的表面高程，可以用來派生出坡度、坡向和流域屬性。光譜值可在衛星影像和航空攝影中表示光反射係數和顏色。

單元值可正可負，可以是整型也可以是浮點型。整數值適合表示類別（離散）數據；浮點值則適合表示連續表面。有關離散和連續數據的附加資訊，請參閱離散和連續數據。在單元中，還可以使用 NoData 值來表示數據缺失。有關 NoData 的資訊，請參閱柵格數據集中的 NoData。

柵格數據的像元值

柵格會被儲存爲有序的像元值列表，例如：80、74、62、45、45、34 等。

各像元所表示區域（或表面）的高和寬都相等，而且在柵格表示的整個表面上佔據相等的部分。例如，表示高程的一個柵格（即，數位高程模型）可能會覆蓋 100 平方千米的區域。如果該柵格中有 100 個像元，則每個像元都將表示等高等寬的 1 平方千米（即，1 km x 1 km）。

像元的大小決定了柵格中要素呈現的粗細程度。

像元的尺寸可大可小，具體可根據柵格數據集所描述的表面，以及表面中要素的表達需要來確定。它可以是平方千米、平方英尺，甚至是平方釐米。像元的大小決定着柵格中圖案或要素呈現的粗細程度。像元大小越小，則柵格將越平滑或越詳細。但是像元數量越多，所需的處理時間會越長，佔據的儲存空間也越大。如果像元大小過大，則可能會出現資訊丟失或精細的圖樣變得模糊的情況。例如，如果像元大小超過道路的寬度，則柵格數據集中便不存在該道路。下圖顯示如何使用不同像元大小的柵格數據集來表示簡單的面要素。

將數據儲存爲柵格具有以下優點：
1、數據結構更加簡單
2、格式更加強大，可進行高階的空間和統計分析
3、可以表示連續表面以及執行表面分析
4、點、線、面和表面都可同樣儲存
5、對複雜數據集也可執行快速疊置

將數據儲存爲柵格具有以下侷限性：例如：
1、像元尺寸具有侷限性可能帶來空間誤差。
2、柵格數據集會非常大佔用更多的磁碟空間，拖慢處理速度。
3、數據重建會損失一定的精度。

柵格數據基本詞彙

您在 ArcGIS 中使用柵格數據時會遇到的一些術語進行概述。這些術語並非全部特定於 ArcGIS。

詞彙	描述
柵格與影像	柵格和影像是兩個經常互相指代的術語。影像是二維的影象表示。它不依賴於波長或遙感裝置，如衛星、航空攝像機或地形感測器。影像可以顯示在螢幕上，也可以列印出來。您可以檢視影像。柵格是描述影像儲存方式的數據模型。柵格可定義組成影像的畫素數（像元數）（以行和列的形式表示）、波段數以及位深度。當您檢視柵格時，您檢視的是該柵格數據的影像。您還可能聽說柵格被稱爲基於像元的數據集，但在 ArcGIS 文件中通常不這樣使用。
像元和畫素	畫素通常會作爲像元的同義詞使用。像元和畫素都是指柵格數據中的最小資訊單位。畫素是影象元素的簡稱，通常用於描述影像，而像元則通常用於描述柵格數據。像元和畫素具有尺寸和值。它們可以表示溫度、土壤型別、高程和真實世界要素（例如，公園、湖泊和建築物）等資訊。
解析度、比例和像元大小	解析度、比例和像元（畫素）大小都是指柵格數據中要素的大小，但又並非如此簡單。例如，有四種類型的解析度：光譜解析度 - 描述用於建立影像的電磁光譜內的波長。時間解析度 - 指對地球表面的同一地點捕獲影像的頻率。輻射解析度 - 描述感測器在電磁光譜的同一部分中對所檢視物件的分辨能力。空間解析度 - 唯一一個與比例和像元大小相關的解析度。空間解析度（也稱爲像元大小）是單個像元所表示的地面上覆蓋的區域的尺寸。比例是地圖（或影像）上的距離或面積與地面上相對應的距離或面積之間的比率或關係，通常以分數或比率的形式表示。空間解析度或像元大小會影響影像在任意比例下所表示的詳細程度。
波段	柵格可包含一個或多個波段。多波段柵格通常稱爲多光譜影象，而具有高達數百個波段的柵格通常稱爲高光譜影象。單波段柵格數據集表示單一現象（如高程）或只表示電磁光譜的一個波長範圍（如黑白航空像片）。波段通常與光譜解析度相關。
柵格格式與柵格型別	柵格格式用於定義畫素的儲存方式，例如，行數和列數、波段數、實際畫素值，以及其他柵格格式特定的參數。柵格型別用於與柵格格式一起標識元數據，例如地理配準、採集日期和感測器型別。
柵格產品	柵格產品的作用是使從特定感測器或數據提供商向您的地圖新增影像變得更加容易，因爲每一柵格產品均有獨特的增強功能的集合以及波段組合以優化您數據的檢視。柵格產品將在目錄中顯示，以代替與特定供應商產品相關的元數據檔案。它是用於生成柵格產品（比如 Landsat 7 或 QuickBird 等衛星影像）的元數據檔案中的資訊。柵格產品會以自身獨有的圖示顯示在目錄中：柵格產品圖示。
渲染	柵格數據集可在地圖中以多種不同的方式進行顯示或渲染。渲染是一個顯示數據的過程。柵格數據集的渲染方式取決於它所包含的數據的型別以及您要顯示的內容。某些柵格包含一個 ArcMap 將自動用於顯示柵格數據的預定義配色方案，即色彩對映表。對於未包含預定義配色方案的柵格，ArcMap 將選擇一種合適的顯示方法，您可根據需要對其進行調整。
功能分類	函數允許您（或軟體）定義將應用於一個或多個柵格的處理功能，但此處理功能不會永久地應用於柵格；而是在存取柵格時動態應用。

	儲存方法：數據模型 瞭解有關柵格數據模型的資訊
柵格數據集	柵格數據集是組織成一個或多個波段的任何有效的柵格格式。每個波段由一系列畫素（單元）陣列組成，每個畫素都有一個值。柵格數據集至少有一個波段。可將多個柵格數據集在空間上附加（鑲嵌）在一起形成一個更大的連續柵格數據集。柵格數據集用 柵格數據集圖標表示：
鑲嵌數據集	鑲嵌數據集是一組以目錄形式儲存並以單個鑲嵌影像或單個影像（柵格）的方式顯示或存取的柵格數據集（影像）。這些集合的總檔案大小和柵格數據集數量都會非常大。新增柵格數據時會根據其柵格型別進行，該型別與柵格格式一起用於標識元數據，例如地理配準、採集日期和感測器型別。鑲嵌數據集中的柵格數據集可以採用本機格式保留在磁碟上，也可在需要時載入到地理數據庫中。可通過柵格記錄以及屬性表中的屬性來管理元數據。通過將元數據儲存爲屬性，可以更方便地管理諸如感測器方向數據等參數，同時也可以提高對選擇內容的查詢速度。鑲嵌數據集用鑲嵌數據集圖標表示：
柵格目錄	柵格目錄是以表格式定義的柵格數據集的集合，其中每個記錄表示目錄中的一個柵格數據集。柵格目錄可以大到包含數千個影像。柵格目錄通常用於顯示相鄰、完全重疊或部分重疊的柵格數據集，而無需將它們鑲嵌爲一個較大的柵格數據集。柵格目錄用 柵格目錄圖標表示：
託管與非託管	在地理數據庫中儲存柵格數據的方式有兩種 - 託管或非託管。託管柵格數據集儲存在地理數據庫中，而非託管柵格數據集則不是。

在 ArcGIS 中使用鑲嵌
術語「鑲嵌」在 ArcGIS 中的使用方式有很多種：
鑲嵌數據集 - 地理數據庫中的數據模型，用於管理以目錄形式儲存並以鑲嵌影像方式檢視的柵格數據集（影像）的集合（上文進行了詳細介紹）。這種鑲嵌數據集使用鑲嵌數據集工具集中的地理處理工具進行建立和編輯。
鑲嵌的影像 - 在檢視中顯示的來自兩個或多個影像的影像或柵格。通常，在您檢視鑲嵌數據集時會看到
鑲嵌 - 描述將多個柵格數據集組合或合併在一起的操作。有關詳細資訊，請參閱什麼是鑲嵌。
鑲嵌柵格數據集或柵格數據集（鑲嵌）表示同一含義。二者均爲柵格數據集。詞語「鑲嵌」或「鑲嵌的」用於描述柵格數據集的建立情況，是已建立還是將要建立。

什麼是GDAL？

GDAL(Geospatial Data Abstraction Library)是一個用於柵格數據操作的庫，是開源地理空間基金會（Open Source Geospatial Foundation，OSGeo）的一個專案。

GDAL是一個操作各種柵格地理數據格式的庫。包括讀取、寫入、轉換、處理各種柵格數據格式（有些特定的格式對一些操作如寫入等不支援）。它使用了一個單一的抽象數據模型就支援了大多數的柵格數據（GIS對柵格，向量，3D數據模型的抽象能力實在令人歎服）。

當然除了柵格操作， GDAL從2.0起整合了另一個用於向量數據操作的庫（OGR）；這個庫還同時包括 了操作向量數據的另一個有名的庫ogr（ogr這個庫另外介紹），這樣這個庫就同時具備了操作柵格和向量數據的能力。

fiona一個用於讀寫向量空間數據檔案的Python包，是由Sean Gillies等人寫的，其底層是C++開發的OGR庫（一個開源GIS庫）。

但是ogr這個庫，只是簡單的向量數據的讀寫操作，如若需要一些複雜的幾何操作，那麼就需要配合GEOS庫了，很多開源的GIS軟體（如QGIS）是利用GEOS庫開發的。

GDAL/OGR庫是用C/C++寫的，很多GIS產品都使用了GDAL/OGR庫，包括ESRI的ArcGIS、Google Earth、GRASS GIS等。

GDAL提供了多種程式語言的API，包括Python、Perl、C#以及Java等，因此，在Python中呼叫GDAL的API函數時，底層執行的是C/C++編譯的二進制檔案。

GDAL是基於C/C++開發的，直接使用pip安裝有時會有問題，可從http://www.lfd.uci.edu/~gohlke/pythonlibs/網站上下載GDAL的wheel檔案進行安裝。

python環境如何安裝gdal參考：

一個涉及命令列操作的教學，可以使用附帶的gdal的命令工具

GDAL/OGR Python文件（https://gdal.org/python/）

GDAL/OGR操作手冊

如何對柵格數據進行讀取

GDAL的Open(filename)函數用於讀柵格數據，函數返回Dataset物件。
GDAL目前支援約100種格式的柵格數據讀取，包括ERDAS Imagine、ENVI、GRASS、GeoTIFF、HDF4、HDF5、TIFF、JPEG、JPEG2000、PNG、GIF、BMP等。詳細資訊可檢視https://gdal.org/drivers/raster/index.html。

GDAL支援的柵格數據檔案

通過Dataset物件，可以瞭解柵格數據的基本資訊，如行列數、波段數、座標轉換參數等。
通過Dataset物件可以返回每個波段數據（Band物件），進一步瞭解每個波段數據的資訊。
Dataset物件和Band物件都可以轉換成陣列，通常情況下，柵格數據是基於陣列進行操作。

獲取柵格數據基本資訊
Dataset物件的RasterYSize、RasterXSize和RasterCount屬性分別返回柵格數據的行數、列數和波段數。

from osgeo import gdal
ds = gdal.Open("D:/20200309/sanwei/kuangshan/landsat8001.img")
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount
print(f"rows:{rows}")
print(f"cols:{cols}")
print(f"bands:{bands}")

Dataset物件的GetGeoTransform()方法返回柵格數據的座標轉換參數，即行列座標與空間座標的轉換參數。
返回的值是個元組，共有6個元素，其中，第一和第四個元素爲左上角像元的x和y座標，第二和第六個元素爲x和y方向的比例尺，第三和第五元素爲x和y方向旋轉角度。

from osgeo import gdal
ds = gdal.Open("D:/20200309/sanwei/kuangshan/landsat8001.img")
ds.GetGeoTransform()

Dataset物件的GetProjection()方法返回柵格數據的空間參照系統資訊（WKT文字）。

from osgeo import gdal
ds = gdal.Open("D:/20200309/sanwei/kuangshan/landsat8001.img")
ds.GetProjection()

** 返回波段數據**

柵格數據通常有多個波段， Dataset物件的GetRasterBand(index)方法將返回某個波段的數據物件（Band物件），如ds.GetRasterBand(1)返回第一個波段的數據物件（注意：第一波段是1，不是0）。

通過Band物件可以返回波段數據的相關資訊：
GetMinimum()和GetMaximum() 方法返回波段數據的最小值和最大值。
ComputeStatistics()方法返回波段數據的統計資訊，包括最小值、最大值、平均值和標準偏差，approx_ok參數表示是否抽樣統計，值爲False表示統計所有柵格。

柵格數據轉陣列

Dataset物件和Band物件都有ReadAsArray()方法用於把柵格數據轉換成陣列。
如果是單波段柵格數據，ReadAsArray()函數返回的是二維陣列（rows,columns）；如果是多波段柵格數據，ReadAsArray()函數返回的是三維陣列（bands,rows,columns）。

from osgeo import gdal
ds = gdal.Open("D:/20200309/sanwei/kuangshan/landsat8001.img")
ds_array = ds.ReadAsArray()
print(ds_array.shape)
band1 = ds.GetRasterBand(1)
band1_array = band1.ReadAsArray()
print(band1_array.shape)

對單波段柵格數據，ReadAsArray()函數返回的二維陣列可直接利用imshow()函數進行顯示。

對多波段的柵格數據，由於ReadAsArray()函數返回的三維陣列是bands、rows、columns順序，而imshow()函數則要求三維陣列的順序是rows、columns、bands，因此需要對傳入的數據進行轉置，即transpose(1,2,0)，轉置後，原先的第一維變成第三維，原先的第二維變成第一維，原先的第三維變成第二維。

此外，對三維陣列，imshow()函數要求輸入數據是0～1的浮點數或0～255整數。

import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal
ds = gdal.Open( "D:/20200309/sanwei/kuangshan/landsat8001.img" )
band1_array = ds.GetRasterBand(1).ReadAsArray()
mean = np.mean(band1_array)
std = np.std(band1_array)
plt.imshow(band1_array,
           vmin=mean-std*2,
           vmax=mean+std*2,
           cmap=plt.cm.gray)

//對單波段的柵格數據進行拉伸顯示

如果柵格數據的值範圍很大，但同時大多數柵格的值又集中在一個很小範圍，顯示影象就很難區分柵格之間的差異，通常需要對影象的拉伸顯示，一種簡單的拉伸方法是以柵格數據的平均值減去2倍標準差爲最小值，以柵格數據的平均值加上2倍標準差爲最大值。

import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal
ds = gdal.Open( "D:/20200309/sanwei/kuangshan/landsat8001.img" )
array = ds.ReadAsArray()
array_543 = array[[4,3,2]].transpose(1,2,0)
mean = np.mean(array_543)
std = np.std(array_543)
vmin=mean-std*2
vmax=mean+std*2
clip = np.clip(array_543,vmin,vmax)
a = np.array((clip-vmin)/(vmax-vmin)*255,dtype=int)
plt.imshow(a)
//三波段數據顯示

在轉換成陣列時，可以設定轉換數據的範圍，關鍵字參數xoff和yoff爲起始柵格的列和行的值，xsize和ysize（Band物件是win_xsize和win_ysize）爲轉換數據的列數和行數，預設情況下是轉換整個數據。

柵格數據行列號和地理座標相互轉換

** 行列座標與地理座標相互轉換**

利用座標轉換參數可以進行行列座標與空間座標的相互轉換。

行列座標轉空間座標（柵格中心點空間座標）的計算公式如下：
x = columnpixelWidth + originX – pixelWidth/2
y = linepixelHeight + originY – pixelHeight/2

line和column爲行列座標（行列號）；pixelWidth和pixelHeigt爲x方向比例尺和y方向比例尺；originX和originY爲左上角的x和y座標

from osgeo import gdal
def Pixel2world(geotransform, line, column):
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    x = column*pixelWidth + originX - pixelWidth/2
    y = line*pixelHeight + originY - pixelHeight/2
    return(x,y)
ds = gdal.Open( "D:/20200309/sanwei/kuangshan/landsat8001.img" )
print(ds.GetProjection())
geotransform = ds.GetGeoTransform()
line = 300;column = 200
x,y = Pixel2world(geotransform, line, column)
print(f"x:{x}")
print(f"y:{y}")

空間座標轉行列座標：

空間座標轉行列座標（柵格中心點空間座標）的計算公式如下：

column = int((x – originX) / pixelWidth) + 1
line = int((y – originY) / pixelHeight) + 1

from osgeo import gdal
def world2Pixel(geotransform, x, y):
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    line = int((y-originY)/pixelHeight)+1
    column = int((x-originX)/pixelWidth)+1
    return (line,column)
ds = gdal.Open( "D:/20200309/sanwei/kuangshan/landsat8001.img" )
geotransform = ds.GetGeoTransform()
x = 359880.0;y = 3468780.0
line,column = world2Pixel(geotransform, x, y)
print(f"line:{line}")
print(f"column:{column}")
/*根據空間座標計算行列座標*/

這裏的空間座標 可以理解成 行列號

如何寫入到柵格數據檔案

GDAL寫柵格數據檔案的方式是先產生一個空白（值均爲0）的柵格數據檔案，然後分波段把陣列數據寫入到柵格數據。

柵格數據寫入到檔案的步驟如下：
利用gdal的GetDriverByName(driver_name)函數返回Driver物件（和檔案型別相關）。
利用Driver物件建立一個空的柵格數據檔案，同時返回Dataset物件。
從Dataset物件中返回Band物件，利用Band物件的WriteArray(array)方法把陣列寫入到柵格數據的對應波段中。

利用Dataset物件的FlushCache()方法把記憶體中的數據寫入到檔案中。寫入檔案後，一般需要把Dataset物件設定爲None，否則檔案會被鎖定。

GetDriverByName(driver_name)函數中的driver_name是檔案型別的名稱，如TIFF/GeoTIFF檔案的名稱爲GTiff，Erdas Imagine檔案的名稱爲HFA，ENVI檔案的名稱爲ENVI。
詳細的柵格檔名稱可參考https://gdal.org/drivers/raster/index.html。

GetDriverByName(driver_name)函數中的driver_name是檔案型別的名稱，如TIFF/GeoTIFF檔案的名稱爲GTiff，Erdas Imagine檔案的名稱爲HFA，ENVI檔案的名稱爲ENVI。
詳細的柵格檔名稱可參考https://gdal.org/drivers/raster/index.html。

數據型別決定了柵格值的範圍，如數據型別爲GDT_Byte，則柵格值的範圍是0～255；如要儲存的數據柵格值的範圍是0～65535，則數據型別應該是GDT_UInt16。

%matplotlib inline
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy as np
driver = gdal.GetDriverByName("GTiff")
outdata = driver.Create("c:/tmp/new_image.tif",
                        xsize=300,ysize=200,bands=3)
outdata.FlushCache()
outdata = None
img = plt.imread("c:/tmp/new_image.tif")
plt.imshow(img)

%matplotlib inline
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
driver = gdal.GetDriverByName("GTiff")
outdata = driver.Create("c:/tmp/new_image.tif",
                        ysize=200,xsize=300,bands=3)
for i in range(3):
    outband=outdata.GetRasterBand(i+1)  
    in_array = np.random.rand(200,300)*255
    outband.WriteArray(in_array.astype(int))
outdata.FlushCache()
outdata = None
img = plt.imread("c:/tmp/new_image.tif")
plt.imshow(img)

/*對多個波段數據進行回圈賦值*/

讓我一起來看一個官網的的向量裁剪柵格的例子：

下面 下麪的修改後的指令碼考慮到了這一點，併爲裁剪的Geotiff設定了正確的x，y偏移量。 注意，在以下範例中，我們假設您已安裝Python Imaging Library。

from osgeo import gdal, gdalnumeric, ogr, osr
import Image, ImageDraw
import os, sys
gdal.UseExceptions()


# This function will convert the rasterized clipper shapefile
# to a mask for use within GDAL.
def imageToArray(i):
    """
    Converts a Python Imaging Library array to a
    gdalnumeric image.
    """
    a=gdalnumeric.fromstring(i.tostring(),'b')
    a.shape=i.im.size[1], i.im.size[0]
    return a

def arrayToImage(a):
    """
    Converts a gdalnumeric array to a
    Python Imaging Library Image.
    """
    i=Image.fromstring('L',(a.shape[1],a.shape[0]),
            (a.astype('b')).tostring())
    return i

def world2Pixel(geoMatrix, x, y):
  """
  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
  the pixel location of a geospatial coordinate
  """
  ulX = geoMatrix[0]
  ulY = geoMatrix[3]
  xDist = geoMatrix[1]
  yDist = geoMatrix[5]
  rtnX = geoMatrix[2]
  rtnY = geoMatrix[4]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / xDist)
  return (pixel, line)

#
#  EDIT: this is basically an overloaded
#  version of the gdal_array.OpenArray passing in xoff, yoff explicitly
#  so we can pass these params off to CopyDatasetInfo
#
def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):
    ds = gdal.Open( gdalnumeric.GetArrayFilename(array) )

    if ds is not None and prototype_ds is not None:
        if type(prototype_ds).__name__ == 'str':
            prototype_ds = gdal.Open( prototype_ds )
        if prototype_ds is not None:
            gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
    return ds

def histogram(a, bins=range(0,256)):
  """
  Histogram function for multi-dimensional array.
  a = array
  bins = range of numbers to match
  """
  fa = a.flat
  n = gdalnumeric.searchsorted(gdalnumeric.sort(fa), bins)
  n = gdalnumeric.concatenate([n, [len(fa)]])
  hist = n[1:]-n[:-1]
  return hist

def stretch(a):
  """
  Performs a histogram stretch on a gdalnumeric array image.
  """
  hist = histogram(a)
  im = arrayToImage(a)
  lut = []
  for b in range(0, len(hist), 256):
    # step size
    step = reduce(operator.add, hist[b:b+256]) / 255
    # create equalization lookup table
    n = 0
    for i in range(256):
      lut.append(n / step)
      n = n + hist[i+b]
  im = im.point(lut)
  return imageToArray(im)

def main( shapefile_path, raster_path ):
    # Load the source data as a gdalnumeric array
    srcArray = gdalnumeric.LoadFile(raster_path)

    # Also load as a gdal image to get geotransform
    # (world file) info
    srcImage = gdal.Open(raster_path)
    geoTrans = srcImage.GetGeoTransform()

    # Create an OGR layer from a boundary shapefile
    shapef = ogr.Open(shapefile_path)
    lyr = shapef.GetLayer( os.path.split( os.path.splitext( shapefile_path )[0] )[1] )
    poly = lyr.GetNextFeature()

    # Convert the layer extent to image pixel coordinates
    minX, maxX, minY, maxY = lyr.GetExtent()
    ulX, ulY = world2Pixel(geoTrans, minX, maxY)
    lrX, lrY = world2Pixel(geoTrans, maxX, minY)

    # Calculate the pixel size of the new image
    pxWidth = int(lrX - ulX)
    pxHeight = int(lrY - ulY)

    clip = srcArray[:, ulY:lrY, ulX:lrX]

    #
    # EDIT: create pixel offset to pass to new image Projection info
    #
    xoffset =  ulX
    yoffset =  ulY
    print "Xoffset, Yoffset = ( %f, %f )" % ( xoffset, yoffset )

    # Create a new geomatrix for the image
    geoTrans = list(geoTrans)
    geoTrans[0] = minX
    geoTrans[3] = maxY

    # Map points to pixels for drawing the
    # boundary on a blank 8-bit,
    # black and white, mask image.
    points = []
    pixels = []
    geom = poly.GetGeometryRef()
    pts = geom.GetGeometryRef(0)
    for p in range(pts.GetPointCount()):
      points.append((pts.GetX(p), pts.GetY(p)))
    for p in points:
      pixels.append(world2Pixel(geoTrans, p[0], p[1]))
    rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
    rasterize = ImageDraw.Draw(rasterPoly)
    rasterize.polygon(pixels, 0)
    mask = imageToArray(rasterPoly)

    # Clip the image using the mask
    clip = gdalnumeric.choose(mask, \
        (clip, 0)).astype(gdalnumeric.uint8)

    # This image has 3 bands so we stretch each one to make them
    # visually brighter
    for i in range(3):
      clip[i,:,:] = stretch(clip[i,:,:])

    # Save new tiff
    #
    #  EDIT: instead of SaveArray, let's break all the
    #  SaveArray steps out more explicity so
    #  we can overwrite the offset of the destination
    #  raster
    #
    ### the old way using SaveArray
    #
    # gdalnumeric.SaveArray(clip, "OUTPUT.tif", format="GTiff", prototype=raster_path)
    #
    ###
    #
    gtiffDriver = gdal.GetDriverByName( 'GTiff' )
    if gtiffDriver is None:
        raise ValueError("Can't find GeoTiff Driver")
    gtiffDriver.CreateCopy( "OUTPUT.tif",
        OpenArray( clip, prototype_ds=raster_path, xoff=xoffset, yoff=yoffset )
    )

    # Save as an 8-bit jpeg for an easy, quick preview
    clip = clip.astype(gdalnumeric.uint8)
    gdalnumeric.SaveArray(clip, "OUTPUT.jpg", format="JPEG")

    gdal.ErrorReset()


if __name__ == '__main__':

    #
    # example run : $ python clip.py /<full-path>/<shapefile-name>.shp /<full-path>/<raster-name>.tif
    #
    if len( sys.argv ) < 2:
        print "[ ERROR ] you must two args. 1) the full shapefile path and 2) the full raster path"
        sys.exit( 1 )

    main( sys.argv[1], sys.argv[2] )

看到這裏，相信小夥伴基本上都會操作自己的數據了吧；

注：本文主要目的是對自己學習知識的總結和分享，以便於一些GIS入門的小夥伴參考； 其中也有自己學習的時候：借鑑的圖片和素材，如有侵權小編會及時更正，還望諒解； 本文如若有錯誤，還請在文末流言，小編會及時更正，懇請批評指正。

柵格數據介紹部分，參考：什麼是柵格數據
https://desktop.arcgis.com/zh-cn/arcmap/10.3/manage-data/raster-and-images/what-is-raster-data.htm