Python地学分析 — GDAL通过矢量裁剪遥感图像

欢迎关注博主的微信公众号：“智能遥感”。

该公众号将为您奉上Python地学分析、爬虫、数据分析、Web开发、机器学习、深度学习等热门源代码。

本人的GitHub代码资料主页（持续更新中，多给Star，多Fork）：

https://github.com/xbr2017

CSDN也在同步更新：

https://blog.csdn.net/XBR_2014

“ 通过矢量图层裁剪遥感栅格数据，就可以得到我们的感兴趣区，需要用到掩膜、矢量栅格化等知识。”

今天的遥感之美封面图—国家公园雅乌河, 热带雨林亚马逊。那里有各种各样的生物资源,也有让人流连忘返的热带雨林风光,还有各种值得研究的自然现象。不如咱们现在就通过遥感图像一睹为快吧！公园的名字来自于河流中一种名为"雅乌"的鱼类,这种鱼个头肥大,而且肉质鲜美,是当地人酷爱的佳肴。而且公园内最大的一条河也叫"雅乌",可见当地人是多么热爱这种鱼！

此图由Landsat7的ETM+传感器获得（2008年5月4日）。巴西的黑水河是亚马逊河最大的支流，河流穿过该国的雅乌国家公园（JNP）。当雨季来临时，水位上涨，河道中可见的岛屿会被河水淹没。

雅乌国家公园位于巴西亚马逊地区，南纬1º00-3º00到西经61º30-64º00直接的地域。最初它是雅乌河和黑水河的汇合处，后来，雅乌河的右岸逐渐拓展，一直延伸到卡拉宾那尼河流入口处。再后来，公园又继续顺着卡拉宾那尼河的右岸向上扩展，最后再一次与黑河相遇。公园的海拔高度介于0至200米之间。公园中的森林与连续分布于辽阔的亚马孙中部平原上的森林相连，四周是由地势较低的黑河环绕，从而形成这里一道独特的风景景观。由于保护了亚马逊森林以及雅乌河流域盆地大量有代表性的典型物种，为研究很多自然现象对该地区生物物种多样性的影响提供了独特的机会，于2000年被列入世界遗产名录。

矢量裁剪栅格

在上一节里面，我们讲了遥感图像的拼接。大家有没有想过，如果遥感图像拼接的空间范围比我们感兴趣的研究区域稍微大一些，或者我们的研究区域是一个不规则的形状，如北京市、江苏省等区域，那该怎么办？这时候，我们就需要感兴趣区域的矢量图来裁剪对应该区域的遥感图像。图像裁剪是一个形象地比喻，也是遥感数据预处理常用的手段之一，好比将一个模具放在一块布上，然后按照这个模具的边缘进行裁剪，最终得到与模具一样形状的布，遥感图像裁剪也是如此。

对遥感图像进行裁剪，需要以下主要步骤：

（1）将遥感图像作为数值载入；

（2）读取shp矢量文件；

（3）将shp矢量文件转化为带有地理参考信息的栅格图（矢量转栅格）；

（4）将栅格化之后的shp文件作为“模具”，舍弃模具以外的遥感图像像元；

（5）使用一步中的掩膜文件对遥感图像进行裁剪；

（6）删除不在边框范围内的像元；

（7）最后将裁剪后的数据保存为Geotiff文件输出。

这里需要解释一下什么叫做掩膜，它的英文单词是Mask，对，没错！就叫马赛克。如果连马赛克都不知道的话，那你只能去请教苍老师了！！！可能是翻译老师考虑到不让学生们浮想联翩吧，所以就叫掩膜啦。下面用一幅示意图来看看图像掩膜的原理：

将原始图像与掩膜图像相乘之后，得到掩膜后的效果图

掩膜前的遥感图像（上面叠加了一个面状矢量图层）：

矢量裁剪遥感图像代码实现：

# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2019/1/2 18:06'import operator
from osgeo import gdal, gdal_array, osr
import shapefiletry:import Imageimport ImageDraw
except:from PIL import Image, ImageDraw# 用于裁剪的栅格数据
raster = r'D:\FalseColor.tif'
# 用于裁剪的多边形shp文件
shp = r'D:\hancock'
# 裁剪后的栅格数据
output = r'D:\clip'def image2Array(i):"""将一个Python图像库的数组转换为一个gdal_array图片"""a = gdal_array.numpy.fromstring(i.tobytes(), 'b')a.shape = i.im.size[1], i.im.size[0]return adef world2Pixel(geoMatrix, x, y):"""使用GDAL库的geomatrix对象((gdal.GetGeoTransform()))计算地理坐标的像素位置"""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) / abs(yDist))return (pixel, line)# 将数据源作为gdal_array载入
srcArray = gdal_array.LoadFile(raster)
# 同时载入gdal库的图片从而获取geotransform
srcImage = gdal.Open(raster)
geoTrans = srcImage.GetGeoTransform()
# 使用PyShp库打开shp文件
r = shapefile.Reader("{}.shp".format(shp))
# 将图层扩展转换为图片像素坐标
minX, minY, maxX, maxY = r.bbox
ulX, ulY = world2Pixel(geoTrans, minX, maxY)
lrX, lrY = world2Pixel(geoTrans, maxX, minY)
# 计算新图片的尺寸
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)
clip = srcArray[:, ulY:lrY, ulX:lrX]
# 为图片创建一个新的geomatrix对象以便附加地理参照数据
geoTrans = list(geoTrans)
geoTrans[0] = minX
geoTrans[3] = maxY
# 在一个空白的8字节黑白掩膜图片上把点映射为像元绘制市县
# 边界线
pixels = []
for p in r.shape(0).points:pixels.append(world2Pixel(geoTrans, p[0], p[1]))
rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
# 使用PIL创建一个空白图片用于绘制多边形
rasterize = ImageDraw.Draw(rasterPoly)
rasterize.polygon(pixels, 0)
# 使用PIL图片转换为Numpy掩膜数组
mask = image2Array(rasterPoly)
# 根据掩膜图层对图像进行裁剪
clip = gdal_array.numpy.choose(mask, (clip, 0)).astype(gdal_array.numpy.uint8)
print clip.max()
# 将NDVI保存为tiff文件
gdal_array.SaveArray(clip, "{}.tif".format(output),format="GTiff", prototype=raster)

裁剪后的遥感图像：

Python地学分析 — GDAL通过矢量裁剪遥感图像相关推荐

Python地学分析 — GDAL将多个遥感图像叠加保存为tif文件
欢迎关注博主的微信公众号:"智能遥感". 该公众号将为您奉上Python地学分析.爬虫.数据分析.Web开发.机器学习.深度学习等热门源代码. 本人的GitHub代码资料主页(持续 ...
Python地学分析 — GDAL对遥感影像重投影
欢迎关注博主的微信公众号:"智能遥感". 该公众号将为您奉上Python地学分析.爬虫.数据分析.Web开发.机器学习.深度学习等热门源代码. 本人的GitHub代码资料主页(持续 ...
Python地学分析 — 矢量数据集介绍
欢迎关注博主的微信公众号:"智能遥感". 该公众号将为您奉上Python地学分析.爬虫.数据分析.Web开发.机器学习.深度学习等热门源代码. 什么叫矢量数据,在这里,小编先给大家 ...
Python地学分析 — 建立矢量数据缓冲区 06
欢迎关注博主的微信公众号:"智能遥感". 该公众号将为您奉上Python地学分析.爬虫.数据分析.Web开发.机器学习.深度学习等热门源代码. Python的小伙伴们,你们好!上一 ...
Python地学分析 — 地理空间参考系介绍
欢迎关注博主的微信公众号:"智能遥感". 该公众号将为您奉上Python地学分析.爬虫.数据分析.Web开发.机器学习.深度学习等热门源代码. 本人的GitHub代码资料主页(持续 ...
ENVI裁剪遥感图像问题
想把上海市裁剪出来,请问大家这个上海市的面状数据导入之后为什么跟遥感图像的边界相差这么多?这是正常的吗?是面状数据不对还是软件操作的问题呢?
Python遥感图像处理应用篇(十九)：GDAL +numpy批量对遥感图像外围背景值进行处理
1.问题描述最近下载了一些遥感影像数据,这些数据都包含大量的外围背景数据,如下图所示: 外围背景值都为0值. 本文描述的是采用python批量处理外围背景,只保留最小外围背景区域. 如下图: 2.实 ...
自主开发的遥感图像数据处理系统
功能包括: 遥感图像的控制点数据库管理,入库.检索.查找等基于航空影像的三维地形重建,生成DEM.DOM数据可见光.SAR.红外图像,相互之间的配准.融合序列图像(可见光.红外)的配准镶嵌基本 ...
Python使用GDAL矢量裁剪栅格，设置背景值为空白（已解决）
一.使用gdal.Warp gdalwarp 实用程序是一种图像拼接.重投影和扭曲实用程序.该程序可以重新投影到任何支持的投影.如果图像是带有控制信息的"原始"图像,也可以存储原始 ...
Python遥感图像处理应用篇(二十二)：Python+GDAL 批量等距离裁剪影像-续
之前写过一篇按照指定行列号数量来进行影像等距离裁剪的博客,链接如下: Python遥感图像处理应用篇(二十二):Python+GDAL 批量等距离裁剪影像_空中旋转篮球的博客-CSDN博客_pytho ...

Python地学分析 — GDAL通过矢量裁剪遥感图像

Python地学分析 — GDAL通过矢量裁剪遥感图像相关推荐

最新文章

热门文章