Option to create a vector field map from saved maps (OperatoR)

[arunge]

• add option to create a vector field map from saved maps using: https://r-graphics.org/recipe-miscgraph-vectorfield. Standard display options for generated map then apply.

[arunge]

@isomemo I have a question regarding the logic that should be used to define the vector field.

In the example from the link there is data that contains speed components in each direction:

library(gcookbook) # For the isabel data set
isabel
#>                x        y      z        vx        vy          vz         t
#> 1      -83.00000 41.70000  0.035        NA        NA          NA        NA
#> 2      -83.00000 41.55571  0.035        NA        NA          NA        NA
#> 3      -83.00000 41.41142  0.035        NA        NA          NA        NA
#> 156248 -62.12625 24.09679 18.035 -11.39709 -5.315139 0.009657148 -66.99567
#> 156249 -62.12625 23.95251 18.035 -11.37965 -5.275015 0.040921956 -67.00032
#> 156250 -62.12625 23.80822 18.035 -12.16637 -5.435891 0.030216325 -66.98057
#>           speed
#> 1            NA
#> 2            NA
#> 3            NA
#>  ...<156,244 more rows>...
#> 156248 12.57555
#> 156249 12.54281
#> 156250 13.32552

The function geom_segment requires values for xend and yend:

islice <- filter(isabel, z == min(z))

ggplot(islice, aes(x = x, y = y)) +
       geom_segment(aes(xend = x + vx/50, yend = y + vy/50),
                    size = 0.25)   # Make the line segments 0.25 mm thick

• To create a map we currently do plot latitude, longitude, estimates .
• If we have two maps and an operation, e.g. +,
  • we first create the resulting matrix for the selected operation that contains latitude, longitude, estimateAfterApplyingOperation
  • and than we do plot the data to the map.

What values should be used for xend and yend here if we want to apply geom_segment? We can also discuss this in our next meeting.

[isomemo]

@arunge we should discuss to better understand the problem!

[arunge]

@isomemo As discussed, you want to check if there is a different package that already contains the calculation to determine the length of the arrows (they should describe e.g. spread over time). So we will return to this issue later.

[isomemo]

@arunge one way of doing this is to use a D8 algorithm that for each cell checks which of its immediate 8 neighbours has the biggest difference (positive or negative depending on user selection). The arrow in the cell will point in the direction the selected neighbour with a magnitude corresponding the difference in value between the two cells and the distance (square root of the sum of the two pixel sizes). I could not find an R implementation for this? Perhaps someone at INWT can help? Otherwise I will write up the code.

[arunge]

@isomemo Regarding the vector field map in OperateR:

• I will check how (where) to extract the content of a cell

The coloured maps are combining

data1 <- data.frame(
  long1 = c(...),
  lat1 = c(...),
  z1 = c(...)
)
data2 <- data.frame(
  long2 = c(...),
  lat2 = c(...),
  z2 = c(...)
)

by finding the closest mean of z2 for coordinates (lat1, long1) and then doing the operation on z1 & closest_mean and keeping coordinates (lat1, long1)

• calculation of the closest mean is here:
  

  DSSM/R/01-plotMap.R

  Lines 2745 to 2759 in 890412b

  	closestMean <- sapply(1:nrow(XPred1), function(x)
  	XPred2$Est[which.min((XPred1$Longitude[x] - XPred2$Longitude) ^ 2 +
  	(XPred1$Latitude[x] - XPred2$Latitude) ^ 2)])
  	closestSd <- sapply(1:nrow(XPred1), function(x)
  	XPred2$Sd[which.min((XPred1$Longitude[x] - XPred2$Longitude) ^ 2 +
  	(XPred1$Latitude[x] - XPred2$Latitude) ^ 2)])
  	closestDist <- sapply(1:nrow(XPred1), function(x)
  	min((XPred1$Longitude[x] - XPred2$Longitude) ^ 2 +
  	(XPred1$Latitude[x] - XPred2$Latitude) ^ 2
  	))
  	closestMean[closestDist > sqrt(2)] <- NA
  	closestSd[closestDist > sqrt(2)] <- NA
  	XPred <- XPred1

• calculation of different operations:
  

  DSSM/R/01-plotMap.R

  Lines 2760 to 2808 in 890412b

  	if(operation == "-"){
  	XPred$Est <- XPred1$Est - closestMean
  	XPred$Sd <- sqrt(XPred1$Sd ^ 2 + closestSd ^ 2)
  	}
  	if(operation == "+"){
  	XPred$Est <- XPred1$Est + closestMean
  	XPred$Sd <- sqrt(XPred1$Sd ^ 2 + closestSd ^ 2)
  	}
  	if(operation == "*"){
  	XPred$Est <- XPred1$Est * closestMean
  	XPred$Sd <- sqrt(XPred1$Sd ^ 2 * closestSd ^ 2 +
  	XPred1$Sd ^ 2 * closestMean ^ 2 +
  	closestSd ^ 2 * XPred1$Est ^ 2)
  	}
  	if(operation == "/"){
  	XPred$Est <- XPred1$Est / closestMean
  	XPred$Sd <- sqrt(((XPred1$Est ^ 2) / (closestMean ^ 2)) *
  	((XPred1$Sd ^ 2) / (XPred1$Est ^ 2) + (closestMean ^ 2) / (closestSd ^ 2)))
  	}
  	if(operation == "mean"){
  	XPred$Est <- (XPred1$Est + closestMean) / 2
  	XPred$Sd <- sqrt(XPred1$Sd ^ 2 + closestSd ^ 2) / sqrt(2)
  	}
  	if(operation == "weightedMean"){
  	VAR <- 1 / (1 / (XPred1$Sd ^ 2) + 1 / (closestSd ^ 2))
  	XPred$Est <- VAR * (XPred1$Est / (XPred1$Sd ^ 2) + closestMean / (closestSd ^ 2))
  	XPred$Sd <- sqrt(VAR)
  	}
  	if(operation == "weight"){
  	constant <- mean(closestMean)
  	closestMean <- closestMean / constant
  	closestSd <- closestSd / constant
  	XPred$Est <- XPred1$Est * closestMean
  	XPred$Sd <- sqrt(XPred1$Sd ^ 2 * closestSd ^ 2 +
  	XPred1$Sd ^ 2 * closestMean ^ 2 +
  	closestSd ^ 2 * XPred1$Est ^ 2)
  	}
  	if(operation == "pMin"){
  	XPred$Est <- pmin(XPred1$Est, closestMean, na.rm = T)
  	XPred1$Est[is.na(XPred1$Est)] <- -Inf
  	XPred$Sd <- sapply(1:length(XPred1$Sd), function(l) c(XPred1$Sd[l], closestSd[l])[which.min(c(XPred1$Est[l], closestMean[l]))])
  	XPred$Sd[is.na(XPred$Est)] <- NA
  	}
  	if(operation == "pMax"){
  	XPred$Est <- pmax(XPred1$Est, closestMean, na.rm = T)
  	XPred1$Est[is.na(XPred1$Est)] <- Inf
  	XPred$Sd <- sapply(1:length(XPred1$Sd), function(l) c(XPred1$Sd[l], closestSd[l])[which.max(c(XPred1$Est[l], closestMean[l]))])
  	XPred$Sd[is.na(XPred$Est)] <- NA
  	}

• the result is returned here:
  

  DSSM/R/01-plotMap.R

  Lines 2810 to 2812 in 890412b

  	XPred$IntLower <- XPred$Est - 1.96 * XPred$Sd
  	XPred$IntUpper <- XPred$Est + 1.96 * XPred$Sd
  	return(XPred)

[arunge]

@isomemo Just to write down where we left today:

Above I summarized how operateR is calculating the combination of two maps.

We will need a different logic to create vector field maps.

For plotting we can use e.g. geom_segment(), a ggplot and sf (https://r-spatial.org/r/2018/10/25/ggplot2-sf-2.html)

ggplot(data, aes(x = x, y = y)) +
geom_segment(aes(xend = xend, yend = yend)

That is, we need values for xend and yend.

You mentioned a method that one can use for calculating vector fields. Can you give me the name of the method you just suggested?