mlt / QGIS-Processing-tools

Custom QGIS processing scripts
http://mlt.github.io/QGIS-Processing-tools
2 stars 0 forks source link

Channel migration tool not working in QGIS 3 and/or with long centerlines #1

Closed alvitello closed 5 years ago

alvitello commented 5 years ago

Hi, I am getting the following message when trying to run Channel Migration:

Performing alignment...
Error in dtw(with(from.df, data.frame(x, y, Curvature * Curvature_multiplier)), :
No warping path exists that is allowed by costraints
Calls: lapply -> FUN -> dtw
Execution halted

It seem like it depends on the centerline shapefile I am using: some of them work just fine, while some others don't. However, I cannot figure out what the differences are, and so why the code doesn't work sometimes.

Polyline features (see attached file) are single-part and oriented in the same way. I am using QGIS Wien 2.8.9 (didn't get the tools to work on more recent versions).

Best

CentLine.zip

mlt commented 5 years ago

Just thought to get back to you ASAP so you can get started… back in the days I put constraints on window/envelope to minimize computation time. There is some bad (?) heuristics. Apparently it was looking for migrations within very narrow vicinity (1189/10 points) and it may be problematic at times as in this case where maximum shift is 137.

> range( with(alignment, index2-index1) )
[1]   0 137
> nrow(from.df)
[1] 1189

You can try to comment out these two lines. I tend to think I either should skip that as such or use something proportional to log(nrow). Let me know how it goes.

Also note that you can use it from R directly without QGIS

# install.packages(c('sp', 'rgdal'))

library(sp)
library(rgdal)

Original <- readOGR('t0_line.shp')
Final <- readOGR('t1_line.shp')
Years <- 1
spar <- .4
Curvature_multiplier <- 1e5
Step_pattern <- 'symmetricP05'
Step <- 25

source('rscripts/Channel migration.rsx')

writeOGR(Migration, '.', "migration", driver='ESRI Shapefile')

I'll see into newer QGIS compatibility.

mlt commented 5 years ago

It works okay with QGIS 3.6.3-Noosa on Windows 10 for me as is. You ought to download and install Processing R Provider plugin. Note that by default it wants to use Austrian (at ?) R mirror and local folder with R packages instead of global. You might want to check those options. Also script files should go into %APPDATA%\QGIS\QGIS3\profiles\default\processing\rscripts instead of ~/.qgis2/processing, if you are on Windows.

alvitello commented 5 years ago

Thank you, I've got it working in QGIS 3.6.3.-Noosa on Win7 by ignoring the two lines you suggested (without that, I still get the same error message) and by changing the script file folder. Thanks a lot, this is really an excellent tool!

mlt commented 5 years ago

Nice! I am glad to hear it is working for you. Just a heads up that this is purely geometrical approach as you might have seen. There is no accounting whatsoever for any underlying physical processes. It won't work when there are extensive meandering changes like oxbow lake formation. There was no proper review and/or testing but given its simplicity I'm not sure it is necessary. Said that, I would really appreciate if you could drop a line here or by e-mail on how you ended up using this in your work and/or research so I could use that in success story/testimony if you do not mind. Thanks!

ara92mc commented 5 years ago

Hello :)! I also have issues, but I am using R/Rstudio. And what I get is a little line (the "migration" file attached). This is my code:

work_dir <- "C:/Users/gmallma/Documents/Hmcenterline"
setwd(work_dir)
library(sp)
library(rgdal)
library(dtw)

Original <- readOGR('C:/Users/gmallma/Documents/Hmcenterline/h_1987_mc.shp')
Final <- readOGR('C:/Users/gmallma/Documents/Hmcenterline/h_1993_mc.shp')
Years <- 1
spar <- .4
Curvature_multiplier <- 1e5
Step_pattern <- 'symmetricP05'
Step <- 25
source("C:/Users/gmallma/Documents/Hmcenterline/Channel migration.rsx")
writeOGR(Migration, '.', "migration", driver='ESRI Shapefile')

I've also tried to do as suggested before by commenting on those 2 code lines, and it didn´t work. So I don´t really know what I'm doing wrong.

rivers.zip

mlt commented 5 years ago

There are few things off with your data. First and foremost, you don't want to use geographical coordinates. You need to project your layers probably to UTM Zone 18S to get meaningful numbers like meters per year for erosion rate. Also, probably, you would want Years <- 1993-1987 (6). Note that it is a very long river. You will likely run out of memory trying to allocate things with Step <- 25. Either increase Step for more coarse step between migration lines. Or split manually into smaller segments.

...
Original0 <- readOGR('h_1987_mc.shp')
Final0 <- readOGR('h_1993_mc.shp')

Original <- spTransform(Original0, CRS("+proj=utm +zone=18 +South +datum=WGS84"))
Final <- spTransform(Final0, CRS("+proj=utm +zone=18 +South +datum=WGS84"))

Years <- 6
Step <- 100
...
writeOGR(Migration, '.', "migration", driver='ESRI Shapefile', overwrite_layer=TRUE)

I hope this will get you going.

ara92mc commented 5 years ago

Thank you so much! I was just going to do that (the UTM files). It worked fine!

mlt commented 5 years ago

I highly suggest to split very long centerlines into shorter segments. This way you would be able to use small Step values to get better overall accuracy and stats per stream bend without getting out of memory errors. For this, you can guesstimate relatively straight and/or unchanged parts of a river and place some points into a new layer that you can use later to segment original and final centerlines. Then use the following code to cut centerlines at vertices close to points you chose.

split_points0 <- readOGR('split_points.shp')
split_points <- spTransform(split_points0, CRS("+proj=utm +zone=18 +South +datum=WGS84"))

library(rgeos)
split_long_stream <- function(stream, points) {
  pt.dist <- sort( gProject(stream, points) )
  coords <- coordinates(stream)[[1]][[1]]

  pt.idx <- 1
  cumdist <- 0
  ll <- list()
  vertex1 <- 1
  for(j in seq_len(nrow(coords))[-1]) {
    cumdist <- cumdist + as.numeric( dist(coords[(j-1):j,]) )
    if (cumdist >= pt.dist[pt.idx]) {
      ll <- c(ll, Lines( Line(coords[vertex1:j,]), pt.idx ))
      vertex1 <- j
      pt.idx <- pt.idx + 1
      if (pt.idx > length(pt.dist)) break
    }
  }
  ll <- c(ll, Lines( Line(coords[(j+1):nrow(coords),]), pt.idx ))
  spl <- SpatialLines(ll, proj4string=CRS(proj4string(stream)))
  SpatialLinesDataFrame(spl, data.frame(id=1:pt.idx))
}

spl.df <- split_long_stream(Original, split_points)
writeOGR(spl.df, '.', "Original_split", driver='ESRI Shapefile', overwrite_layer=TRUE)

spl.df <- split_long_stream(Final, split_points)
writeOGR(spl.df, '.', "Final_split", driver='ESRI Shapefile', overwrite_layer=TRUE)

Perhaps the entire process can be automated with the first coarse step used to locate good candidate locations for splitting, followed by some weighted clustering, then the second run with the smaller step on pre-segmented centerlines. I'm not sure how to quantify the criteria for splitting yet.

ara92mc commented 5 years ago

Hi, so I have a new error, actually the code is working pretty fine, but I don´t know why do I get this results. Heres, this is my code: work_dir <- "C:/Users/gmallma/Documents/Hmcenterline" setwd(work_dir) library(sp) library(rgdal) library(dtw)

Original0 <- readOGR('C:/Users/gmallma/Documents/Hmcenterline/h_2017_mc.shp') Final0 <- readOGR('C:/Users/gmallma/Documents/Hmcenterline/h_2011_mc(1).shp')

Original <- spTransform(Original0, CRS("+proj=utm +zone=18 +South +datum=WGS84")) Final <- spTransform(Final0, CRS("+proj=utm +zone=18 +South +datum=WGS84"))

Years <- 6 spar <- 0.4 Curvature_multiplier <- 1e5 Step_pattern <- 'symmetricP05' Step <- 500

source("C:/Users/gmallma/Documents/Hmcenterline/Channel migration.rsx")

writeOGR(Migration, '.', "Migra_2011_2017", driver='ESRI Shapefile', overwrite_layer = TRUE)

_____

With the modifications you suggested :) I will attached the results I got, as well as the code.

migracion.zip