前面一篇文章()讲述了通过 geopandas 库实现对子区域数据的分类统计,说白了也就是如何根据一个 shp 数据对另一个 shp 数据进行切割。本篇作为上一篇内容的姊妹篇讲述如何采用优雅的方式根据一个 shp 数据对一个栅格影像数据进行切割。废话不多说,直接进入主题。
本方案涉及以下技术点:
本文基本涉及以上技术。另,最近 Github 貌似被墙了,所以你懂的。推荐使用 Lantern,请自行百度之。
为什么叫优雅的切割,其实我这里倒不是卖弄文字,主要是为了与 Gdal 的方式相区别。传统的方式可以采用 Gdal 命令行进行一点点的手动处理,稍微智能化一点可以在 python 程序中发送控制台语句的方式调用 gdal 命令。作为程序员我们都是想采用最简单、最不需要手工操作、看上去最舒服的方式。所以我这里称其为优雅的方式。
我们大致需要经历读取影像、投影转换、读取 shp、切割、显示等几个步骤。下面逐一介绍。
采用 rasterio 进行影像读取。代码如下:
- import rasterio as rio
- band = rio.open(path)
非常简单,只要传入影像的路径即可。rasterio 支持 tif、hdf 格式(亲测)。
这是比较难的一块,需要注意很多细节。代码如下:
- from rasterio.warp import (reproject,RESAMPLING, transform_bounds,calculate_default_transform as calcdt)
- affine, width, height = calcdt(src.crs, dst_crs, src.width, src.height, *src.bounds)
- kwargs = src.meta.copy()
- kwargs.update({
- 'crs': dst_crs,
- 'transform': affine,
- 'affine': affine,
- 'width': width,
- 'height': height,
- 'geotransform':(0,1,0,0,0,-1) ,
- 'driver': 'GTiff'
- })
- dst = rio.open(newtiffname, 'w', **kwargs)
- for i in range(1, src.count + 1):
- reproject(
- source = rio.band(src, i),
- destination = rio.band(dst, i),
- src_transform = src.affine,
- src_crs = src.crs,
- dst_transform = affine,
- dst_crs = dst_crs,
- dst_nodata = src.nodata,
- resampling = RESAMPLING.bilinear)
首先使用 calculate_default_transform 函数计算投影转换后的影像尺寸以及仿射变换参数。其中 src 表示原始影像,src.crs 可以取到原始投影,dst_crs 位定义的目标投影,与上一篇中介绍 shp 投影变换时基本一致。
然后计算投影后的 tiff 元数据信息。src.meta.copy() 读出原始元数据信息并进行拷贝,kwargs.update 将原始元数据更新为目标元数据。
dst = rio.open(newtiffname,'w', **kwargs) 打开一个新的影像其模式 w 表示写入。
最后循环原始影像的所有波段,逐一进行投影变换并写入新的影像。其参数一目了然,不再赘述。
上一个影像的整体截图,以与下述切割后的效果进行对比。
这在上一篇文章中也已经做了详细描述,不再赘述,需要强调的时此处也需要将 shp 进行投影转换,使其与我们要处理的影像一致,所以简单的方式就是直接读取影像的投影信息,将 shp 数据转换到此投影,详情请参考。
我们要对一个完整的影像进行切割,可以分为两步。首先将 shp 数据转换为 geojson,然后使用 rasterio 进行切割。
rasterio 进行切割时需要传入的时 geojson 对象,而不是普通的 GeoSeries 对象,所以我们需要进行一步转换。代码如下:
- from geopandas import GeoSeries
- features = [shpdata.geometry.__geo_interface__]
其中 shpdata 为读出的 shp 数据对象。如果我们想要获取 shp 中的某条空间数据而不是全部,可以采用如下方式:
- from geopandas import GeoSeries
- features = [GeoSeries(shpdata.geometry[i]).__geo_interface__]
其中 i 表示的是取出元素的序号,最后都要采用 [] 将结果变成数组,因为 rasterio 最后需要传入一个数组参数。
其实有了前面的准备这一步也就变的简单了,直接调用 rio.mask.mask 函数,该函数返回该栅格数据与 features 相交部分的数组结果以及变换信息。代码如下:
- import rasterio.mask
- out_image, out_transform = rio.mask.mask(src, features, crop=True, nodata=src.nodata)
- out_meta = src.meta.copy()
- out_meta.update({"driver": "GTiff",
- "height": out_image.shape[1],
- "width": out_image.shape[2],
- "transform": out_transform})
- band_mask = rasterio.open(newtiffname, "w", **out_meta)
- band_mask.write(out_image)
其中 src 为切割前的影像数据,features 为上一步得到的 shp 数据转换后的 geojson,crop 表示是否对原始影像进行切割,如果为 True 表示将该 geojson 的外界框以外的数据全部删除,既缩小原始影像的大小,只保留外界框以内部分,nodata 表示无值数据,凡是 geojson 外部的数据都会转换成此值。
后面的基本与投影转换后的一致,根据切割的结果生成一个新的影像数据。这样我们就实现了根据 shp 数据对遥感影像进行切割。效果如下:
本文所介绍的技术可以用于对全国的影像数据进行分省切割,或者省的影像数据进行县市切割等。同理与上一篇文章一致的是凡是这种处理子区域的方式都可以采用此技术。当然本文没有介绍如何对遥感影像进行处理,其实非常简单,当我们读出影像数据之后,其就是一个 numpy 的 array 对象,已经变成了纯数学问题,处理完之后只需要附加投影等信息写入新的 tiff 文件即可。
来源: http://www.cnblogs.com/shoufengwei/p/6437934.html