Python+GDAL栅格数据基本操作

什么是栅格数据？

最简形式的栅格由按行和列（或格网）组织的像元（或像素）矩阵组成，其中的每个像元都包含一个信息值（例如温度）。栅格可以是数字航空像片、卫星影像、数字图片或甚至扫描的地图。

为何将数据存储为栅格？

在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