[R-sig-Geo] raster directions to vectorial lines (Stream to Feature)

[Mauricio Zambrano-Bigiarini]

On 22/03/13 05:31, Gavan McGrath wrote:
> Hi Mauricio,
>
> Below is a simple pair of functions for plotting a flow network over an image which should achieve what you were after.  You may want to play around with how your flow direction numbers 1 -9 match up with those in the lookupfidr function. Its a bit cluncky but gets the job done.
>
> The function fdirPlot takes as input a matrix imData the background image, and a matrix fdirs of flow direction integers.
>
> Good luck.

Thank you very much Gavan for you valuable help, and apologise for my 
late reply.

With slight modifications to the code you gave me, I managed to do what 
I want on the screen device:

------------- START reproducible example ----------------
lookupfdir <-function(fdir) {
    switch(paste(fdir,"",sep=""),
        '1' = c(-1,-1),
        '4' = c(-1,0),
        '7' = c(-1,1),
        '8' = c(0,1),
        '9' = c(1,1),
        '6' = c(1,0),
        '3' = c(1,-1),
        '2' = c(0,-1),
        '5' = c(0,0),
            c(-2,-2)
    )
} # 'lookupfdir' end

fdirPlot <- function(imData,fdirs) {
   dimDat<- dim(imData)
   Ly    <- dimDat[1] - 1
   Lx    <- dimDat[2] - 1
   centres.x <- seq(0,1, length.out=dimDat[2])
   centres.y <- seq(0,1, length.out=dimDat[1])
   image(t(imData),axes=FALSE)
   for (i in centres.x) { # x direction
     for (j in centres.y) { # y direction
     dxdy <- lookupfdir(fdirs[i+1,j+1])
     dx   <- dxdy[1]/Lx
     dy   <- dxdy[2]/Ly
     if (dxdy[1] != -2)
      lines( rbind( c(i, j), c(i + dx, j + dy) ) )
     } # IF end
   } # FOR j end
} # 'fdirPlot' end

# Testing
imData <- matrix(rnorm(30), ncol=3, nrow=2)
fdirs  <- matrix(9, ncol=3, nrow=2)

fdirPlot(imData, fdirs)
------------- END reproducible example ----------------


However, I would like to do the same as before but over a raster file:


# creating the two rasters from the scratch
library(raster)
imData         <- raster(nrows=2, ncols=3, xmn=0, xmx=10,  ymn=0, ymx=5)
values(imData) <- rnorm(6)

fdirs     <- raster(nrows=2, ncols=3, xmn=0, xmx=10,  ymn=0, ymx=5)
values(r) <- 9

However, I do not know how to plot the lines over the raster and then 
save the resulting lines as a shapefile.

Do you -or somebody else- have any advice about how to do it ?

Thank in advance for your help,

Mauricio

-- 
=================================================
Water Resources Unit
Institute for Environment and Sustainability (IES)
Joint Research Centre (JRC), European Commission
TP 261, Via Enrico Fermi 2749, 21027 Ispra (VA), IT
webinfo    : http://floods.jrc.ec.europa.eu/
=================================================
DISCLAIMER:
"The views expressed are purely those of the writer
and may not in any circumstances be regarded as sta-
ting an official position of the European Commission"
=================================================
Linux user #454569 -- Ubuntu user #17469
=================================================
"It is not enough to have knowledge;
one must also apply it" (Goethe)
>
>
> lookupfdir <-function(fdir) {
>     switch(paste(fdir,"",sep=""),
>         '1' = c(-1,-1),
>         '2' = c(-1,0),
>         '3' = c(-1,1),
>         '4' = c(0,1),
>         '5' = c(1,1),
>         '6' = c(1,0),
>         '7' = c(1,-1),
>         '8' = c(0,-1),
>         '9' = c(0,0),
>             c(-2,-2)
>     )
> }
>
> fdirPlot <- function(imData,fdirs) {
>    dimDat<- dim(imData)
>    image(imData,axes=FALSE)
>    for (i in 1:dimDat[1]) {
>    for (j in 1:dimDat[2]) {
>    dxdy <- lookupfdir(fdirs[i,j])
>    if (dxdy[1] != -2) {
>     lines(rbind(c(i-0.5,j-0.5)/dimDat[1],c(i-dxdy[1]-0.5,j-dxdy[2]-0.5)/dimDat[1]))
>     } else {
>            if (i==1) {
>              lines(rbind(c(i-0.5,j-0.5)/dimDat[1],c(i-1.5,j-0.5)/dimDat[1]))
>            } else {
>                  if (i==dimDat[1]) {
>                      lines(rbind(c(i-0.5,j-0.5)/dimDat[1],c(i+1-0.5,j-0.5)/dimDat[1]))
>                   } else {
>                          if (j==1) {
>                            lines(rbind(c(i-0.5,j-0.5)/dimDat[1],c(i-0.5,j-1-0.5)/dimDat[1]))
>                          } else {
>                                 if (j==dimDat[2]) {
> 		                 lines(rbind(c(i-0.5,j-0.5)/dimDat[1],c(i-0.5,j+1-0.5)/dimDat[1]))
>                                 }
>                            }
>                       }
>               }
>        }
>    }
>    }
> }
>
> Gavan McGrath
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>