专栏首页点滴积累使用Python以优雅的方式实现根据shp数据对栅格影像进行切割

使用Python以优雅的方式实现根据shp数据对栅格影像进行切割

一、前言

       前面一篇文章(使用Python实现子区域数据分类统计)讲述了通过geopandas库实现对子区域数据的分类统计,说白了也就是如何根据一个shp数据对另一个shp数据进行切割。本篇作为上一篇内容的姊妹篇讲述如何采用优雅的方式根据一个shp数据对一个栅格影像数据进行切割。废话不多说,直接进入主题。

二、涉及到的技术

本方案涉及以下技术点:

  1. geopandas:已经在上一篇文章中简单介绍。
  2. numpy:这是一个开源的数据分析处理库,非常高效、简洁。
  3. rasterio:这是一个开源的影像处理库,非常好用,基本涵盖了常用的功能。也可以配合numpy进行数据计算。
  4. datashader:这是一个开源的大数据可视化库,可以进行遥感影像、矢量数据的可视化。其基于bokeh,bokeh是一个通用的可视化工具,有兴趣的可以参考github,我之前采用Scala语言对其进行了简单的封装,请参考使用bokeh-scala进行数据可视化以及使用bokeh-scala进行数据可视化(2)

本文基本涉及以上技术。另,最近Github貌似被墙了,所以你懂的。推荐使用Lantern,请自行百度之。

三、优雅切割

       为什么叫优雅的切割,其实我这里倒不是卖弄文字,主要是为了与Gdal的方式相区别。传统的方式可以采用Gdal命令行进行一点点的手动处理,稍微智能化一点可以在python程序中发送控制台语句的方式调用gdal命令。作为程序员我们都是想采用最简单、最不需要手工操作、看上去最舒服的方式。所以我这里称其为优雅的方式。

       我们大致需要经历读取影像、投影转换、读取shp、切割、显示等几个步骤。下面逐一介绍。

3.1 读取影像

       采用rasterio进行影像读取。代码如下:

import rasterio as rio

band = rio.open(path)

       非常简单,只要传入影像的路径即可。rasterio支持tif、hdf格式(亲测)。

3.2 投影转换

       这是比较难的一块,需要注意很多细节。代码如下:

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表示写入。

       最后循环原始影像的所有波段,逐一进行投影变换并写入新的影像。其参数一目了然,不再赘述。

       上一个影像的整体截图,以与下述切割后的效果进行对比。

3.3 读取shp

       这在上一篇文章中也已经做了详细描述,不再赘述,需要强调的时此处也需要将shp进行投影转换,使其与我们要处理的影像一致,所以简单的方式就是直接读取影像的投影信息,将shp数据转换到此投影,详情请参考使用Python实现子区域数据分类统计

3.4 切割

       我们要对一个完整的影像进行切割,可以分为两步。首先将shp数据转换为geojson,然后使用rasterio进行切割。

3.4.1 shp数据转换为geojson

       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最后需要传入一个数组参数。

3.4.2 使用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文件即可。

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 安装ClouderaManager以及使用ClouderaManager安装分布式集群的若干细节

    目录 前言 整体介绍 分步安装介绍 总结 一、前言        周末干了近四十个小时中间只休息了五个小时终于成功安装了ClouderaManager以及分布式...

    魏守峰
  • linux下快速列出文件列表的方法

    前言 这两天碰到一个很棘手的问题,需要读取出ubuntu系统中某个目录下所有文件,由于服务器中存储的文件实在太多,导致此过程效率十分低下,动辄需要等待一个小时之...

    魏守峰
  • geotrellis使用(二十五)将Geotrellis移植到spark2.0

    目录 前言 升级spark到2.0 将geotrellis最新版部署到spark2.0(CDH) 总结 一、前言        事情总是变化这么快,前面刚写了一...

    魏守峰
  • 文娱大数据:这只猴子今年特别火!

    <数据猿导读> 已近年尾,2015年的文娱大数据也随之新鲜出炉,回首这一年,我国电影行业可谓喜忧参半,多部影片屡屡刷新票房纪录,只可惜有些影片却是口碑欠佳。要问...

    数据猿
  • 目前学术界最先进的数据包调度器介绍!

    随着链路速度的提高和CPU速度缩放速度的降低,软件中的数据包调度会导致较低的精度和较高的CPU利用率。通过将数据包调度卸载到诸如NIC之类的硬件,可以潜在地克服...

    网络交换FPGA
  • 无道工具网二维码生成接口(API)

    想为博客增加手机扫描二维码阅读的功能,但网上的接口不一定哪天就失效了。再加上自己的工具站

    无道
  • 克隆TikTok后,小扎新晋全球第三大富豪,净资产超1000亿美元

    Facebook 的首席执政官扎克伯格的净资产超过1000亿美元,剽窃 TikTok 的小扎竟然走上了「人生巅峰」!

    新智元
  • 安全成SDN发展阻碍 需加快与SDx技术融合对接

    编者按:安全问题是SDN商用化进程停滞不前的原因之一,一边是想要颠覆传统的改革派,一边是谨小慎微的保守派,双方就安全问题展开持续斗争。但是SDN已经是一种必然,...

    SDNLAB
  • 十问大数据安全分析(大数据安全的小船怎样才能不翻?)

    导语:人类的生产生活每天都在产生大量的数据,并且产生的速度越来越快。新的攻击手段层出不穷,需要检测的数据越来越多,现有的分析技术不堪重负。 ? 安全数据的数量...

    灯塔大数据
  • 开源软件许可

    因为日常工作中用到了,一些开源的产品,每个产品说明中,会有一些开源许可的介绍,各种名字,不很理解其中的含义。

    bisal

扫码关注云+社区

领取腾讯云代金券