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

今天开始着手改造pyris 1.0.
文章地址&＃xff1a;https://doi.org/10.1016/J.ENVSOFT.2018.03.028

Monegaglia&＃xff0c;2018

segment

segment_all( landsat_dirs, geodir, config, maskdir, auto_label&＃61;None )>>>segment( landsat, geodir, maskdir, config, clipdir&＃61;None, auto_label&＃61;None )

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

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

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

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

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

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

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

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

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


改进6之后把大气校正加进去

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

改造后&＃xff1a;pyris.segment( lds, geodir, maskdir, prjdir, cf, clipdir&＃61;clipdir, auto_label&＃61;None )

clean_masks

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

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

skeletonize

skeletonize_all( maskdir, skeldir, config )

vectorize

vectorize_all( geodir, maskdir, skeldir, config, axisdir, use_geo&＃61;True )

这是个函数使用一个迭代算法来计算中心线&＃xff1a;

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

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

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

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

migration_rates

migration_rates( axisfiles, migdir, columns&＃61;(0,1), show&＃61;False, pfreq&＃61;1 )

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

bars_detection

bars_detection( landsat_dirs, geodir, axisdir, migdir, bardir, show&＃61;False )

Sample

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

boundary_detection

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

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

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

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

Choose_windows

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

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

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

Make_windows

geofile,maskfile,windowsflie&＃xff0c;size&＃xff0c;flow_from
geotransf&＃xff08;地理转换参数&＃xff09;&＃xff0c;

shp_export

之后改名为axis2shp
把提取axis的中心线输出为线shp和点shp&＃xff08;带河宽、曲率等字段&＃xff09;&＃xff0c;或者是多段线要素。

mask2shp

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

shp2mask

shp2mask(shpdir,maskdir,tempdir)

Mk_Gif

gifpath&＃61;Mk_Gif(clipdir)

misc即miscellaneous缩写&＃xff0c;这个模块里面放杂七杂八的函数

misc.save_Trans

pcsize&＃61;GeoTransf[&＃39;PixelSize&＃39;]
y&＃61;GeoTransf[&＃39;X&＃39;]
x&＃61;GeoTransf[&＃39;Y&＃39;]
lx&＃61;GeoTransf[&＃39;Ly&＃39;]
ly&＃61;GeoTransf[&＃39;Lx&＃39;]


misc.load_Trans

misc.save_Proj

misc.load_Proj

misc.pnts2shpline

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

misc.pnts2shppoly

同上

misc.setProj

定义投影&＃xff0c;待写

misc.projTrans

对shp文件进行重投影&＃xff0c;改变到指定坐标系

misc.read_xml

读取xml文件的坐标

misc.Crs_sample

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

misc.Reach_sample

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

2. 用crs去spilt bd_line&＃xff0c;得到各个子河段的bd_line
   

3. 遍历各个子河段的2条边界&＃xff0c;构建子河段的面

4. 多个子河段的面放在一起得到一个meshgrid&＃xff0c;它是一个multipolygon

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

crs(line1,line2)

求2条线段的交点

rmove_loops(line)

去掉一个多段线的loop部分

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

1. 影像格式的问题
   主要有uint16&＃xff0c;float64
2. 函数适用范围

河型

对象

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

适用范围

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

河流中心线提取

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

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

河流中心线摆动计算

滩体形态分析

河流拓扑结构提取

断面采样

子河段采样

河道边界提取

影像数据集制作

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

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

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

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


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

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

Gee-python的使用

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

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

1. 影像日期的信息从文件夹名得到&＃xff0c;在usgs&＃xff1a;
   
   但国内从gscloud下载的文件命名如下&＃xff1a;
   
2. 不同波段影像的读取

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

le7_off
le7_on、lt5在gscloud中波段名后面会带个0&＃xff08;B10、B20…&＃xff09;

S2&＃xff1a;clean_masks的使用
skeldir的手动去刺&＃xff0c;这样加快下一步vectorization的速度。

S3&＃xff1a;水陆边界计算

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