热门标签 | HotTags
当前位置:  开发笔记 > 编程语言 > 正文

【原创】利用Python进行河流遥感处理的PyRIS软件开发

今天开始着手改造pyris1.0.文章地址:https:doi.org10.1016J.ENVSOFT.2018.03.028Monegaglia,2

沱沱河断面及河段采样
今天开始着手改造pyris 1.0.
文章地址:https://doi.org/10.1016/J.ENVSOFT.2018.03.028

Monegaglia,2018

PyRIS改造

  • 主函数
    • segment
    • clean_masks
    • skeletonize
    • vectorize
    • migration_rates
    • bars_detection
  • 新增主函数
    • Sample
    • boundary_detection
    • Choose_windows
    • Make_windows
    • shp_export
    • mask2shp
    • shp2mask
    • Mk_Gif
  • misc模块
    • misc.save_Trans
    • misc.load_Trans
    • misc.save_Proj
    • misc.load_Proj
    • misc.pnts2shpline
    • misc.pnts2shppoly
    • misc.setProj
    • misc.projTrans
    • misc.read_xml
    • misc.Crs_sample
    • misc.Reach_sample
    • crs(line1,line2)
    • rmove_loops(line)
  • raster模块
  • gee_func模块
  • 一些问题
    • 河型
    • 对象
  • 应用
    • 适用范围
    • 河流中心线提取
    • 河流中心线摆动计算
    • 滩体形态分析
    • 河流拓扑结构提取
    • 断面采样
    • 子河段采样
    • 河道边界提取
    • 影像数据集制作
    • 从谷歌地球获得shp & 坐标系转换
    • numpy数组、栅格、面shp相互转换
    • Gee-python的使用
  • SI


主函数

segment

segment_all( landsat_dirs, geodir, config, maskdir, auto_label=None )>>>segment( landsat, geodir, maskdir, config, clipdir=None, auto_label=None )

选择提取mask的方式,有NDVI(植被),MNDWI(水体),BAR(SWIR矿物),MIX(以NDVI刻画河道边界,求MNDWI和NDVI交集,再和SWIR求并集,最后去噪)。

使用MIX方法专门用来提取河道平面(包括水体和滩体,而不包括江心洲和滩洲),此方法可以拿来提取河流中心线(和拿MNDWI来提取中心线相比,排除了水位的影响),需要注意的是,影像的河道外最好都是植被不然NDVI不能很好地刻画河道边界,如下图:城市面积在河道周围分布很广,这样提取的河道需要调用Clean_masks来根据底图手动删改。
在这里插入图片描述

  1. 提取河道中心线,包括水体和滩体,如果用的影像河道周围是明显植被,可以用MIX方法,但是如果影像河道周围没植被,那就只能退而求其次用MNDWI近似的求中心线.
  2. 提取主槽,只包括水体,需要控制遥感对应水位区间,只用MNDWI来求中心线.

改进1:原作者是用的指定区域裁剪,我改成了用shapefile,shp可以跨景,可以不覆盖景,但需要定义投影坐标系。
这样,你就可以把你所有的Landsat放一个文件夹下面,当你要提取某一区域时,只需要准备好对应的shpfile就可以了。
注:函数用的是大津算法计算阈值,你的陆地部分要足够多,这样形成的双峰才能比较好捕捉合适的阈值,也就是说shp不能只贴着河道范围,多画大一些(3到4倍河道宽),也不要太大因为浪费计算时间。
问题 :虽然可以跨景,很方便,但是当你的shpfile最小外接矩形很大时,裁剪下来的图也会很大,虽然shp之外没有值,但是这样计算量也会很大。是不是很像GEE的处理方法?
解决 :首先得到遥感影像范围;根据生成带投影坐标的面shp;转换面shp的投影坐标系为裁剪shp的投影坐标系;把用来裁剪的shp和范围面shp求交集,得到新shp再用这个shp来裁剪,以减小计算量。

改进2:为该函数添加了每幅图投影信息的输出prjdir

改进3:由于遥感影像命名规则略微不同,原代码读取landsat的函数适用于USGS下载的图,对于gscluod上下载的图会出错,对该部分进行了修改。详见S1。

改进4:对于长江源辫状河流,河道外没有绿色植被,MNDWI对阈值敏感,如果使用全局otsu算法,丰水期主流附近的小心滩会识别不出来,局部otsu算法可以解决这个问题,但局部算法也会有识别不了细小汊道的情况。在阈值计算方法里面除了local,global。

改进5:加入了融合(sharpen),把多光谱波段融合到15m分辨率,再进行后续计算。
sharpen工具用的gdal里面的gdal_pansharpen.py工具,由于scripts文件夹下面没有__init__.py文件,所以得这样使用:

from osgeo.scripts import gdal_pansharpen
gdal_pansharpen.gdal_pansharpen(.....)

在这里插入图片描述
改进6之后把大气校正加进去

最后实现的功能:合并多光谱数据集>L7去条带(L7correction)+融合(pansharpen)>计算指标mask(MNDWI,NDVI,SWIR2)>去除遥感图边缘锯齿状空值>去除细小水体(水田)和水体空洞(船)>给独立水体打上标签>保存mask和geotransf,proj。

改造后:pyris.segment( lds, geodir, maskdir, prjdir, cf, clipdir=clipdir, auto_label=None )

clean_masks

clean_masks( maskdir/skeldir, geodir=None, config=None, file_only=False )
启动一个交互界面,让你可以手动去支流(去刺),如果支流很短,其实可以不用专门删它,因为在后面的vectorize里面它会被当成毛刺spur去掉,就是增加一些运行时间。所以其实这个功能在曲流河基本可以不用,也不会怎么影响计算时间,因为曲流河水体拓扑结构比较简单。mask和skel的去刺见S2。

改进1:给交互界面添加了一个功能,可以鼠标滚轮放大和缩小图像,这样方便精细化的clean mask.
改进2:原型的交互界面底图是B1波段的灰度图,我给改成了RGB真彩色。

skeletonize

skeletonize_all( maskdir, skeldir, config )

vectorize

vectorize_all( geodir, maskdir, skeldir, config, axisdir, use_geo=True )

这是个函数使用一个迭代算法来计算中心线:

  1. 算法根据你给的flow_from参数确定第一个点,然后依次遍历相邻点,每遍历一个点把这个点删去。
  2. 当遍历到节点(junction),开启更深一层迭代,依次遍历与节点连接的点,开启节点循环。例如与节点连接的有三个点,则开启一个次数为三的节点循环,在每一个循环内把其它的两个点删除,该次循环结束把删除的点恢复。
  3. 在节点循环内继续往下遍历点,同时记录这一次节点循环经过的点个数和对应宽度,直到遍历到下一个节点并开启下一层迭代(步骤2),当这一循环内的n重迭代全部完成,记录该汊道的长、宽。
  4. 比较各个节点循环产生的汊道,选择符合条件的那一条汊道:在大于最长汊道长度四分之三(去毛刺)的汊道中最宽的那一条。

应用下来,这个方法似乎不适合辫状河、游荡河段。
改进:这个迭代算法会重复计算一些节点,可以改进

Output:axis是8×N的数组[[x地理坐标],[y地理坐标],[里程s],[水流方向theta],[曲率Cs],[河半宽B],[x像素坐标],[y像素坐标]]。
这个河半宽B是按照skeleton里面各个栅格点的值确定的,指那个点到非水体点的最近距离。
曲率Cs有三种计算方法:
水流方向theta的单位是rad,而且是顺时针,所以当使用math.三角函数时要把theta添上一个负号。

注1:flow_from这个参数的选择非常关键,如果选择错误导致起始点不在两端而在河流的中间段,会导致只能提取部分中心线。

migration_rates

migration_rates( axisfiles, migdir, columns=(0,1), show=False, pfreq=1 )

返回的mig和axis尺寸相同,包括各中心线上每个点的迁移方向dx,dy和迁移强度dz,而且把中心线划分为了几个编号的弯道。

bars_detection

bars_detection( landsat_dirs, geodir, axisdir, migdir, bardir, show=False )

指标说明
dxTemporal sequence (list) of x component of the centerline migration vectors
dyTemporal sequence (list) of y component of the centerline migration vectors
dzTemporal sequence (list) of the magnitude of the centerline migration vectors
CssTemporal sequence (list) of the smoothed planform curvature
BITemporal sequence (list) of the bed indexes
B12Temporal sequence (list) of the next bed indexes (where the current bend goes)
BUDTemporal sequence (list) of the bend upstream-downstream index

新增主函数

Sample

Sample(cleanmask_dir,shpdir,geodir,prjdir,braideddir,config)
可以选择河段采样还是断面采样。

boundary_detection

检测河流的水陆边界,可以是岸线,也可以是滩线。根据flow_from分辨左右。
**法一:**自动化计算方法,用中心线栅格的两个端点二元膨胀来做,最后输出是栅格的边缘线,计算过程见S3。

  1. 缺点:端点膨胀半径比较大的话会比较花时间,因为我用的是两个膨胀半径差一个单位圆求差得到的圆环,可以考虑更好的生成环的方法。而且会有误差,可以设置最小半径
  2. 去刺+挑刺+折点判断+label+取最长2个:需要比较规则的边界

**法二:**使用maskclean类似方法打开一个交互界面,首先计算edgemask,再在界面里面手动删掉上下边界,还没做呢。

问:拓扑结构的link根据端点的联通数平均值来分级,有什么意义?
**法三:**还想了一种自动化计算方法,是检索中心线两个端点一定范围内的点,用叉乘判断是否删除这些点,这个可能比法一快,之后有空试一下。

Choose_windows

输入:二维(如MNDWI,NDVI)或者三维(多波段)的数组,窗口的size,流向flow_from
选择一幅图,根据axisread方法读取河流的axismask,遍历axis上的每一个点,根据size(窗口边长)来掩这个部分的可见光波段,同时打开一个交互窗口,可以根据←、→键移动窗口,如果是想要的窗口,则点空格,程序会根据geotransf来获取窗口四个角点坐标。
输出:输出的是所有的窗口坐标,以npy格式保存。

这些windows可以在gee里面用来获取数据集,用于GANs训练。

想要制作数据集,这里用原作者在segmentation里面的操作,用n个rect来裁剪。首先得要河道中心线,有了中心线,就可以以中心线上的点为中心画一个rect,rects用npy格式保存,文件名要注明是哪个投影坐标系的npy。

Make_windows

geofile,maskfile,windowsflie,size,flow_from
geotransf(地理转换参数),

shp_export

之后改名为axis2shp
把提取axis的中心线输出为线shp和点shp(带河宽、曲率等字段),或者是多段线要素。

mask2shp

mask2shp(maskdir,geodir,prjdir,shpdir,tempdir)

shp2mask

shp2mask(shpdir,maskdir,tempdir)

Mk_Gif

gifpath=Mk_Gif(clipdir)

misc模块

misc即miscellaneous缩写,这个模块里面放杂七杂八的函数

misc.save_Trans

pcsize=GeoTransf['PixelSize']
y=GeoTransf['X']
x=GeoTransf['Y']
lx=GeoTransf['Ly']
ly=GeoTransf['Lx']

misc.load_Trans


misc.save_Proj


misc.load_Proj


misc.pnts2shpline

将多段线的点转为对应坐标系下的shp

misc.pnts2shppoly

同上

misc.setProj

定义投影,待写

misc.projTrans

对shp文件进行重投影,改变到指定坐标系

misc.read_xml

读取xml文件的坐标

misc.Crs_sample

这里我用的shp和axis,基于shapely做的。
之后改:制作Crs时要保留中点的信息,这个中点不一定是2端点中点

misc.Reach_sample


  1. 根据axis生成两条边界线bd_line和断面集crs,生成bd_line时需要判断点在河流左岸还是右岸,用中心点流向向量和中心点垂直射向两边的向量的叉乘来实现判断。我生成bd_line时设置的宽度是2倍河宽,crs的宽度是2.5倍河宽,保证crs可以切到bd_line。

  2. 用crs去spilt bd_line,得到各个子河段的bd_line
    请添加图片描述

  3. 遍历各个子河段的2条边界,构建子河段的面

  4. 多个子河段的面放在一起得到一个meshgrid,它是一个multipolygon

  5. 遍历multipolygon,每个子河段和水体做交集,求水体面积,空洞面积
    同样根据步长对子河段采样,计算每个子河段水域面积、陆地面积、平均河宽


crs(line1,line2)

求2条线段的交点

rmove_loops(line)

去掉一个多段线的loop部分

raster模块
gee_func模块

放一些使用gee会用到的函数

一些问题
  1. 影像格式的问题
    主要有uint16,float64
  2. 函数适用范围

河型


func河型
segmentall
skeletonall
vectorization弯曲,分汊
migration中心线,岸线
bar_detectionall

对象

研究的对象包括河岸线、主流线(中心线)、滩体、江心洲
岸线提取需要影像对应流量够高且不过高而导致漫滩,根据ndvi不靠谱
3. 多数据源适配
不只landsat,还有哨兵之类。

应用

适用范围

研究的对象包括河岸线、主流线(中心线)的摆动。滩体、江心洲的识别。断面,河段采样。
岸线提取需要影像对应流量够高且不过高而导致漫滩,根据ndvi不靠谱

河流中心线提取

用1函数来提取mask,2函数手动删除栅格部分,3函数用来提取中心栅格,3函数提取的中心栅格不满意的话还可以再用2函数删除branch。
4函数根据中心栅格数据,以坐标形式输出河道中心线,使用它的时候确保中心栅格都是你要的河道的中心,因为河道中心栅格可能被大桥或者大坝隔开。

只能提取一个河道的中心线,也就是一条线,不适合处理辫状河道。这条中心线是河道的中心线,所以不考虑分汊。

河流中心线摆动计算


滩体形态分析


河流拓扑结构提取


断面采样


子河段采样


河道边界提取


影像数据集制作


  1. Choose_Windows在本地制作窗口集
  2. 根据窗口集在GEE里面采集数据集
  3. 对于保存在本地的数据集,如果采集对象不止一个影像集(比如同时采集了哨兵2和LC8),那么还有分辨率不同的问题,把它们重采样到自定的分辨率

从谷歌地球获得shp & 坐标系转换

从谷歌地球上画的矢量图可以导出为kml格式,它们的坐标系是WGS84地理坐标系,EPSG代号为4326;
下面把它转换为国家2000坐标系,具体对应的投影坐标系代号要去查,如下图,根据你区域的经度范围选择对应的EPSG代号4523(下面这个表格是我在csdn上下载的)
在这里插入图片描述

from osgeo import osr
#读取坐标
lines = pyris.load(r'C:\Users\23932\Desktop\场区范围.kml')#生成面/线shp
pyris.pnts2shpline(lines,r'D:\DATA\TuoTuo\GIS\New_Shapefile(3).shp',proj=None)
pyris.pnts2shppoly(lines, r'D:\DATA\TuoTuo\GIS\New_Shapefile(4).shp', proj=None)#坐标系转换:转换为目标坐标系
pyris.projTrans( r'D:\DATA\TuoTuo\GIS\New_Shapefile(4).shp', r'D:\DATA\TuoTuo\GIS\Newle1.shp', 32646 , inEPSG = 4326)

numpy数组、栅格、面shp相互转换

统计分析有时候会需要面shp,有时候需要用到数组。
可以全部用gdal实现

Gee-python的使用

ROI获取:需要的是UTM84的经纬度坐标。这个坐标有几个途径可以获得,1)GEE engine里面直接手画roi然后获取坐标。2)在谷歌地球里面画roi,然后读取kml坐标信息。3)arcgis里面根据遥感底图话roi,再读取shp坐标信息并转换为WGS84地理坐标(基本不会用到)。

SI

不影响主函数功能理解的部分函数,改进内容放在这里。
S1: LoadLandsatData函数改造
USGS和gscloud命名不同
pyris读取的内容包括波段影像影像日期

  1. 影像日期的信息从文件夹名得到,在usgs:
    在这里插入图片描述
    但国内从gscloud下载的文件命名如下:
    在这里插入图片描述
  2. 不同波段影像的读取

BGRNIRMIRSWIR
LT5123457
LE7123457
LC8234567

lc8、le7_off在gscloud和usgs上都是一样的内容;

le7_offLE7_off
le7_on、lt5在gscloud中波段名后面会带个0(B10、B20…)
在这里插入图片描述
在这里插入图片描述

S2:clean_masks的使用
skeldir的手动去刺,这样加快下一步vectorization的速度。
在这里插入图片描述
S3:水陆边界计算
在这里插入图片描述
在这里插入图片描述

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

在这里插入图片描述
左右岸没分开是因为分辨率太低,而且河太窄

左右岸没分开是因为这里河宽太窄,分辨率又低,正常的话是分得开的,左岸为1,右岸赋值为2.


推荐阅读
  • iOS超签签名服务器搭建及其优劣势
    本文介绍了搭建iOS超签签名服务器的原因和优势,包括不掉签、用户可以直接安装不需要信任、体验好等。同时也提到了超签的劣势,即一个证书只能安装100个,成本较高。文章还详细介绍了超签的实现原理,包括用户请求服务器安装mobileconfig文件、服务器调用苹果接口添加udid等步骤。最后,还提到了生成mobileconfig文件和导出AppleWorldwideDeveloperRelationsCertificationAuthority证书的方法。 ... [详细]
  • Nginx使用AWStats日志分析的步骤及注意事项
    本文介绍了在Centos7操作系统上使用Nginx和AWStats进行日志分析的步骤和注意事项。通过AWStats可以统计网站的访问量、IP地址、操作系统、浏览器等信息,并提供精确到每月、每日、每小时的数据。在部署AWStats之前需要确认服务器上已经安装了Perl环境,并进行DNS解析。 ... [详细]
  • Linux服务器密码过期策略、登录次数限制、私钥登录等配置方法
    本文介绍了在Linux服务器上进行密码过期策略、登录次数限制、私钥登录等配置的方法。通过修改配置文件中的参数,可以设置密码的有效期、最小间隔时间、最小长度,并在密码过期前进行提示。同时还介绍了如何进行公钥登录和修改默认账户用户名的操作。详细步骤和注意事项可参考本文内容。 ... [详细]
  • VScode格式化文档换行或不换行的设置方法
    本文介绍了在VScode中设置格式化文档换行或不换行的方法,包括使用插件和修改settings.json文件的内容。详细步骤为:找到settings.json文件,将其中的代码替换为指定的代码。 ... [详细]
  • 本文介绍了Java工具类库Hutool,该工具包封装了对文件、流、加密解密、转码、正则、线程、XML等JDK方法的封装,并提供了各种Util工具类。同时,还介绍了Hutool的组件,包括动态代理、布隆过滤、缓存、定时任务等功能。该工具包可以简化Java代码,提高开发效率。 ... [详细]
  • android listview OnItemClickListener失效原因
    最近在做listview时发现OnItemClickListener失效的问题,经过查找发现是因为button的原因。不仅listitem中存在button会影响OnItemClickListener事件的失效,还会导致单击后listview每个item的背景改变,使得item中的所有有关焦点的事件都失效。本文给出了一个范例来说明这种情况,并提供了解决方法。 ... [详细]
  • 本文介绍了C#中生成随机数的三种方法,并分析了其中存在的问题。首先介绍了使用Random类生成随机数的默认方法,但在高并发情况下可能会出现重复的情况。接着通过循环生成了一系列随机数,进一步突显了这个问题。文章指出,随机数生成在任何编程语言中都是必备的功能,但Random类生成的随机数并不可靠。最后,提出了需要寻找其他可靠的随机数生成方法的建议。 ... [详细]
  • 本文介绍了如何使用C#制作Java+Mysql+Tomcat环境安装程序,实现一键式安装。通过将JDK、Mysql、Tomcat三者制作成一个安装包,解决了客户在安装软件时的复杂配置和繁琐问题,便于管理软件版本和系统集成。具体步骤包括配置JDK环境变量和安装Mysql服务,其中使用了MySQL Server 5.5社区版和my.ini文件。安装方法为通过命令行将目录转到mysql的bin目录下,执行mysqld --install MySQL5命令。 ... [详细]
  • 突破MIUI14限制,自定义胶囊图标、大图标样式,支持任意APP
    本文介绍了如何突破MIUI14的限制,实现自定义胶囊图标和大图标样式,并支持任意APP。需要一定的动手能力和主题设计师账号权限或者会主题pojie。详细步骤包括应用包名获取、素材制作和封包获取等。 ... [详细]
  • 本文介绍了在Mac上安装Xamarin并使用Windows上的VS开发iOS app的方法,包括所需的安装环境和软件,以及使用Xamarin.iOS进行开发的步骤。通过这种方法,即使没有Mac或者安装苹果系统,程序员们也能轻松开发iOS app。 ... [详细]
  • 本文讨论了如何使用Web.Config进行自定义配置节的配置转换。作者提到,他将msbuild设置为详细模式,但转换却忽略了带有替换转换的自定义部分的存在。 ... [详细]
  • 本文介绍了闭包的定义和运转机制,重点解释了闭包如何能够接触外部函数的作用域中的变量。通过词法作用域的查找规则,闭包可以访问外部函数的作用域。同时还提到了闭包的作用和影响。 ... [详细]
  • YOLOv7基于自己的数据集从零构建模型完整训练、推理计算超详细教程
    本文介绍了关于人工智能、神经网络和深度学习的知识点,并提供了YOLOv7基于自己的数据集从零构建模型完整训练、推理计算的详细教程。文章还提到了郑州最低生活保障的话题。对于从事目标检测任务的人来说,YOLO是一个熟悉的模型。文章还提到了yolov4和yolov6的相关内容,以及选择模型的优化思路。 ... [详细]
  • 本文介绍了如何清除Eclipse中SVN用户的设置。首先需要查看使用的SVN接口,然后根据接口类型找到相应的目录并删除相关文件。最后使用SVN更新或提交来应用更改。 ... [详细]
  • 本文介绍了使用cacti监控mssql 2005运行资源情况的操作步骤,包括安装必要的工具和驱动,测试mssql的连接,配置监控脚本等。通过php连接mssql来获取SQL 2005性能计算器的值,实现对mssql的监控。详细的操作步骤和代码请参考附件。 ... [详细]
author-avatar
蛋蛋小可爱的诱惑_360
这个家伙很懒,什么也没留下!
PHP1.CN | 中国最专业的PHP中文社区 | DevBox开发工具箱 | json解析格式化 |PHP资讯 | PHP教程 | 数据库技术 | 服务器技术 | 前端开发技术 | PHP框架 | 开发工具 | 在线工具
Copyright © 1998 - 2020 PHP1.CN. All Rights Reserved | 京公网安备 11010802041100号 | 京ICP备19059560号-4 | PHP1.CN 第一PHP社区 版权所有