the path/segment thing

[mdsumner]

Building a table of all box coordinates with each row containing its directed face (.vx0 start, .vx1 end - ordered around the box) is seemingly way too hard from the current arrangement.

It should be much more straightforward. I think the problem is that p0, p1 records the orientation of vertices within a face, but there's no orientation of a face with respect to a box (and the path the linked segments take). This must be implicitl at least in the BGM, so record it explicitly in the tables and make it easy to to do.

noba_file  <- bgmfiles::bgmfiles("Nordic")

u <- "https://raw.githubusercontent.com/cecilieha/NoBA/master/corners_neighbours_nordic.txt"
x <- readr::read_tsv(u)
plot(c(x$Lon1, x$Lon2), c(x$Lat1, x$Lat1), type = "n")
segments(x$Lon1, x$Lat1, x$Lon2, x$Lat2)

library(rbgm)

## BGM is a doubly-connected-edge-list
## here in 'related tables' form
dc_edge_list <- bgmfile(noba_file)

library(dplyr)


box_verts <- dc_edge_list$boxes %>% 
     dplyr::select(.bx0, label) %>% 
     inner_join(dc_edge_list$boxesXverts, ".bx0") %>% 
  inner_join(dc_edge_list$vertices, ".vx0") %>% 
  group_by(.bx0) %>% mutate(irow = row_number()) %>% ungroup()

face_vertex_ind <- match(box_verts$.vx0, dc_edge_list$facesXverts$.vx0)
box_verts$face_id <- dc_edge_list$facesXverts$.fx0[face_vertex_ind]
box_verts$face_order <- dc_edge_list$facesXverts$.p0[face_vertex_ind]


face_box_ind <- match(box_verts$face_id, dc_edge_list$faces$.fx0)
box_verts$Neighbour <- ifelse(dc_edge_list$faces$left[face_box_ind] == box_verts$.bx0,dc_edge_list$faces$right, dc_edge_list$faces$left)

ptcentre <- as.data.frame(boxSpatial(dc_edge_list))
coordinates(ptcentre) <- c("insideX", "insideY")
proj4string(ptcentre) <- CRS(dc_edge_list$extra$projection)
library(mapview)
library(sf)
mapview(faceSpatial(dc_edge_list), lwd = 4, color = "pink")  + 
  mapview(ptcentre)
## note too that from here conversion to Spatial is straightforward
## spbabel::sp(box_verts %>% transmute(object_ = .bx0, branch_ = label, x_ = x, y_ = y, island_ = 1, order_ = irow))


## note that .vx0 refers to the start vertex in terms of start =  x1,y1, end = x2,y2
box_verts %>% transmute(irow, .bx0, x1 = x, y1 = y, .vx0) %>% 
  inner_join(box_verts %>% transmute(irow = irow - 1, x2 = x, y2 = y, .bx0)) %>% distinct(.vx0)


%>% 
  left_join(dc_edge_list$facesXverts %>% dplyr::filter(.p0 == 2))


f## the key entity we want is each face, with four coordinates at each end 
## but grouped by box
face_vertex <- dc_edge_list$facesXverts %>% 
  mutate(.p0 = c("start_vertex", "end_vertex")[.p0]) %>% 
  tidyr::spread(.p0, .vx0) %>% 
  inner_join(dc_edge_list$vertices, c("end_vertex" = ".vx0")) %>% 
  rename(x1 = x, y1 = y) %>%
  inner_join(dc_edge_list$vertices, c("start_vertex" = ".vx0")) %>% 
  rename(x2 = x, y2 = y)
face_vertex[c("Lon1", "Lat1")] <- rgdal::project(as.matrix(face_vertex[c("x1", "y1")]), 
                                                 dc_edge_list$extra$projection, inv = TRUE)
face_vertex[c("Lon2", "Lat2")] <- rgdal::project(as.matrix(face_vertex[c("x2", "y2")]), 
                                                 dc_edge_list$extra$projection, inv = TRUE)
## we have one row per face
plot(c(x$Lon1, x$Lon2), c(x$Lat1, x$Lat1), type = "n")
segments(x$Lon1, x$Lat1, x$Lon2, x$Lat2)
with(face_vertex, segments(Lon1, Lat1, Lon2, Lat2, col = "firebrick", lwd = 2))


## faces is fine but we really want box-focus
face_vertex %>% inner_join(dc_edge_list$facesXboxes, c(".fx0" = "iface"))

[RMORSEcode]

I have the solution most of the way there based on Alex' atlantistools,, the only detail I still cannot get is the 'neighbor', which is wrong in the solution below (which is calculated similarly to your code above). There should be an easy way to merge faces data from 'faces.rm' and the table 'newdf2' but I can't seem to get it to work quite right to yield the proper length of 572.

`
library(atlantistools)
noba_file  <- bgmfiles::bgmfiles("Nordic")
box=load_box(noba_file)

### get details for box vertices, depth and area
  proj_in <- box$projection
  if (!any(grepl(pattern = "[+]", proj_in))) {
    proj_in <- strsplit(proj_in, "[[:space:]]")
    for (i in seq_along(proj_in)) {
      proj_in[[i]] <- paste0("+", proj_in[[i]])
    }
    proj_in <- paste(unlist(proj_in), collapse = " ")
  }
  
  n_boxes <- box$nbox
  n_faces=box$nface
  
  # Extract appearance if boxes.
  result <- list()
  for (i in 1:n_boxes) {
    result[[i]] <- data.frame(box$boxes[[i]]$vert)
    names(result[[i]]) <- c("lat", "long")
    result[[i]]$inside_lat <- box$boxes[[i]]$inside[1]
    result[[i]]$inside_long <- box$boxes[[i]]$inside[2]
    result[[i]]$polygon <- i - 1
    result[[i]]$Area = box$boxes[[i]]$area
    result[[i]]$meanDepth = box$boxes[[i]]$botz
    # result[[i]]$iface=box$boxes[[i]]$iface # face number in list of vertices
    # result[[i]]$ibox=box$boxes[[i]]$ibox # neighbor for same face number
    
  }
  result <- do.call(rbind, result)
  
  # Convert coordinates to map-coordinates!
  lat_long <- proj4::project(result[, 1:2], proj = proj_in, inverse = T)
  result$long <- lat_long$x
  result$lat <- lat_long$y
  lat_long <- proj4::project(result[, 3:4], proj = proj_in, inverse = T)
  result$inside_long <- lat_long$x
  result$inside_lat <- lat_long$y
  


### get details on faces

  r.boxXfaces=list()
  for (i in 1:n_boxes) {
    r.boxXfaces[[i]] <- data.frame(rep(i-1,length(box$boxes[[i]]$ibox)))
    r.boxXfaces[[i]]$iface=box$boxes[[i]]$iface # face number in list of vertices
    r.boxXfaces[[i]]$ibox=box$boxes[[i]]$ibox # neighbor for same face number
  }
  r.boxXfaces <- do.call(rbind, r.boxXfaces)
  colnames(r.boxXfaces)[1]='polygon'
    
r.faces <- list()
  for (i in 1:n_faces) {
    r.faces[[i]]=data.frame(box$faces[[i]][1,'p1'], fix.empty.names = F, row.names = i)
    # names(r.faces[[i]]) <- "Lat1"
    r.faces[[i]]$Lat1 =data.frame(box$faces[[i]][1,'p1'],fix.empty.names = F,row.names = i)
    r.faces[[i]]$Lon1 =data.frame(box$faces[[i]][2,'p1'],fix.empty.names = F,row.names = i)
    r.faces[[i]]$Lat2 =data.frame(box$faces[[i]][1,'p2'],fix.empty.names = F,row.names = i)
    r.faces[[i]]$Lon2 =data.frame(box$faces[[i]][2,'p2'],fix.empty.names = F,row.names =i)
    r.faces[[i]]$length.km =data.frame(box$faces[[i]][1,'length']/1000,fix.empty.names = F,row.names = i)
    r.faces[[i]]$l=data.frame(box$faces[[i]][1,'lr'],fix.empty.names = F, row.names=i)
    r.faces[[i]]$r=data.frame(box$faces[[i]][2,'lr'],fix.empty.names = F, row.names=i)
    
    r.faces[[i]][1]=NULL
    names(r.faces[[i]])=c("Lat1", "Lat2", "Lon1", "Lon2", "length.km", "L", "R")
  }
  result.faces2 <- do.call(rbind, r.faces)
  result.faces=data.frame(matrix(unlist(r.faces), nrow=n_faces, ncol=7, byrow = T))
  result.faces$iface=seq(from=0, to=n_faces-1, by=1)
  
  # Convert coordinates to map-coordinates
  lat_long <- proj4::project(result.faces[, 1:2], proj = proj_in, inverse = T)
  result.faces$Lat1 <- lat_long$y
  result.faces$Lon1 <- lat_long$x
  lat_long <- proj4::project(result.faces[, 3:4], proj = proj_in, inverse = T)
  result.faces$Lat2 <- lat_long$y
  result.faces$Lon2 <- lat_long$x
  result.faces=result.faces[,5:11]
  colnames(result.faces)[1]='length.km'
  colnames(result.faces)[2]='L'
  colnames(result.faces)[3]='R'
  

### Create table similar to Cecilie's
### atlantistools solution from 'result' created above, from modified 'convert_begm_Cecilie'
newdf=result[1,]
newdf[1,]=NA
newdf$Lat2=NA
newdf$Lon2=NA
for (i in 1:length(unique(result$polygon))){
  ii=i-1
  xx=subset(result[result$polygon==ii,])
  for (j in 2:dim(xx)[1]-1){
    xx$Lat2[j]=xx$lat[j+1]
    xx$Lon2[j]=xx$long[j+1]
  }
  xx=xx[-c(dim(xx)[1]),]
  newdf=rbind(newdf, xx)
}
newdf2=newdf[-1,]
colnames(newdf2)
colnames(newdf2)[1]='Lat1'
colnames(newdf2)[2]='Lon1'

### Final details on faces:
faces.rm=merge(r.boxXfaces, result.faces, by='iface', all.x=T)`

## test plot faces by order (black, red, green, blue, cyan, magenta, yellow, gray, (rep),)
tt=newdf2[newdf2$polygon==1,]
plot(c(tt$Lon1, tt$Lon2), c(tt$Lat1, tt$Lat2), type = "n")
for(i in 1:dim(tt)[1]){
segments(tt$Lon1[i], tt$Lat1[i], tt$Lon2[i], tt$Lat2[i], col=i)
}

[RMORSEcode]

This works now, should do the trick!

`

### RM 08172017 script to create 'Neighbors' file in order to parse ROMS data and calculate exchange across box faces
### for Atlantis models based on BGM geometry files

library(atlantistools)
library(bgmfiles)
bgm.file  <- bgmfiles::bgmfiles("Nordic")
box=convert_bgm(bgm.file)

## Modified from 'convert_bgm()' from atlantistools to include necessary extra information
convert_bgm_4ROMS <- function(bgm) {
  box <- load_box(bgm = bgm)
  
  # Get info of projection used! Some models don't use '+' to split their
  # arguments in the projection. So we add them here ;)
  proj_in <- box$projection
  if (!any(grepl(pattern = "[+]", proj_in))) {
    proj_in <- strsplit(proj_in, "[[:space:]]")
    for (i in seq_along(proj_in)) {
      proj_in[[i]] <- paste0("+", proj_in[[i]])
    }
    proj_in <- paste(unlist(proj_in), collapse = " ")
  }
  
  n_boxes <- box$nbox
  n_faces=box$nface
  
  # Extract appearance if boxes.
  result <- list()
  for (i in 1:n_boxes) {
    result[[i]] <- data.frame(box$boxes[[i]]$vert)
    names(result[[i]]) <- c("lat", "long")
    result[[i]]$inside_lat <- box$boxes[[i]]$inside[1]
    result[[i]]$inside_long <- box$boxes[[i]]$inside[2]
    result[[i]]$polygon <- i - 1
    result[[i]]$Area = box$boxes[[i]]$area
    result[[i]]$meanDepth = box$boxes[[i]]$botz
  }
  result <- do.call(rbind, result)
  
  # Convert coordinates to map-coordinates
  lat_long <- proj4::project(result[, 1:2], proj = proj_in, inverse = T)
  result$long <- lat_long$x
  result$lat <- lat_long$y
  lat_long <- proj4::project(result[, 3:4], proj = proj_in, inverse = T)
  result$inside_long <- lat_long$x
  result$inside_lat <- lat_long$y
  return(result)
}


### Take output from convert_bgm_4ROMS, remove duplicate values in each box used to close polygon, rename cols
newdf=result[1,]
newdf[1,]=NA
newdf$Lat2=NA
newdf$Lon2=NA
for (i in 1:length(unique(result$polygon))){
  ii=i-1
  xx=subset(result[result$polygon==ii,])
  for (j in 2:dim(xx)[1]-1){
    xx$Lat2[j]=xx$lat[j+1]
    xx$Lon2[j]=xx$long[j+1]
  }
  xx=xx[-c(dim(xx)[1]),]
  newdf=rbind(newdf, xx)
}
newdf2=newdf[-1,]
colnames(newdf2)
colnames(newdf2)[1]='Lat1'
colnames(newdf2)[2]='Lon1'  

### Find which vertices intersect with which polygons
test1=data.frame(newdf2$Lat1, newdf2$Lon1, newdf2$polygon)
test6=aggregate(test1, by=list(test1[,1], test1[,2]), FUN=unique)

### Match endpoints of polygon segments to polygon vertices
# get polygons for each point of segment, find like connection between points, append neighbor polygon to newdf2
newdf2.Lat1Lon1.in.Group1=match(paste(newdf2$Lat1, newdf2$Lon1, sep=','), paste(test6$Group.1, test6$Group.2, sep=','))
newdf2.Lat2Lon2.in.Group1=match(paste(newdf2$Lat2, newdf2$Lon2, sep=','), paste(test6$Group.1, test6$Group.2, sep=','))
toMatch=data.frame(newdf2.Lat1Lon1.in.Group1, newdf2.Lat2Lon2.in.Group1)
toMatch=data.frame(newdf2.Lat1Lon1.in.Group1, newdf2.Lat2Lon2.in.Group1)
toMatch$poly1=test6$newdf2.polygon[toMatch$newdf2.Lat1Lon1.in.Group1]
toMatch$poly2=test6$newdf2.polygon[toMatch$newdf2.Lat2Lon2.in.Group1]
toMatch$box=newdf2$polygon # original polygon
for (i in 1:dim(toMatch)[1]){
  a=toMatch$poly1[[i]] # boxes for Lat1 Lon1
  b=toMatch$poly2[[i]] # boxes for Lat2 Lon2
  c=Reduce(intersect, list(a,b)) # find common boxes between points
  d=toMatch$box[i] # polygon box number for face segment
  if (length(c) < 2){ # select neighbor box that is not equal to polygon box
    if (c==d){
      e=NA # 
    } else {
      e=c
    }
  }else if (length(c)==2) {
    e=c[which(c!=d)]
  }
  toMatch$Neighbor[i]=e #select neighbor box, not boxface
}
newdf2$Neighbor=toMatch$Neighbor
plot_boxes(data=box) # check to make sure neighbors make sense

[mdsumner]

Cool thanks, hopefully will get a chance to explore soon!