【原创】基于Python的河流之遥感处理（PyRIS）

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

Monegaglia，2018

环境

scikit-image0.19.3
pandas1.5.3

主函数

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这个参数的选择非常关键，如果选择错误导致起始点不在两端而在河流的中间段，会导致只能提取部分中心线。

曲率计算

CurvaturePCS(*args, **kwargs)
======================

Compute the channel centerline's arc-length, inflection angle and curvature
from a PCS. Three methods are available:
  1 - Finite differences
  2 - Guneralp and Rhoads 2007 (requires spline derivatives)
  3 - Schwenk et al. 2015

默认的是1算法

migration_rates

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

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

指标	说明
dx	Temporal sequence (list) of x component of the centerline migration vectors
dy	Temporal sequence (list) of y component of the centerline migration vectors
dz	Temporal sequence (list) of the magnitude of the centerline migration vectors
Css	Temporal sequence (list) of the smoothed planform curvature
BI	Temporal sequence (list) of the bed indexes
B12	Temporal sequence (list) of the next bed indexes (where the current bend goes)
BUD	Temporal sequence (list) of the bend upstream-downstream index

$dz=(d{x}^2+d{y}^2{)}^{0.5}$
Css简化后的曲率.
BUD 的值应该是对应弯道的特征值，按2，-1，0，1,2来的，2表示拐点，0表示定点。

bars_detection

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

新增主函数

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	河型
segment	all
skeleton	all
vectorization	弯曲，分汊
migration	中心线，岸线
bar_detection	all

对象

研究的对象包括河岸线、主流线（中心线）、滩体、江心洲
岸线提取需要影像对应流量够高且不过高而导致漫滩，根据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. 不同波段影像的读取

	B	G	R	NIR	MIR	SWIR
LT5	1	2	3	4	5	7
LE7	1	2	3	4	5	7
LC8	2	3	4	5	6	7

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

le7_off
le7_on、lt5在gscloud中波段名后面会带个0（B10、B20…）

S2：clean_masks的使用
skeldir的手动去刺，这样加快下一步vectorization的速度。

S3：水陆边界计算

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