neuralCity.R

library(rgdal)
library(maptools)
library(spdep)
library(bezier)
library(rgeos)

bati = readOGR("D:/DATAS/APUR", "EMPRISE_BATIE")

coords.total = coordinates(bati)

periodes = unique(bati$C_PERCONST)[order(unique(bati$C_PERCONST))]

makeLineFromCoords <- function(coords, i) {
  Sl1 = Line(coords)
  S1 = Lines(list(Sl1), ID=as.character(i))
  Sl = SpatialLines(list(S1))  
  return(Sl)
}

for (periode in periodes) {
  
  print(periode)
  
  coords = coords.total[bati$C_PERCONST <= periode, ] # WE CUMULATE THE RESULTS
  
  # MINIMUM SPANNING TREE
  nb = tri2nb(coords)
  nbw = nb2listw(nb, zero.policy = FALSE)
  mst.bh <- mstree(nbw,5)
  
  # WE BUILD LINES
  synapses = vector(mode="list")
  for (i in 1:nrow(mst.bh)) {

    from = mst.bh[i, 1]
    to = mst.bh[i, 2]
    
    d = sqrt((coords[from, 1] - coords[to, 1])^2 + (coords[from, 2] - coords[to, 2])^2)
    
    midpoint = c((coords[from, 1] + coords[to, 1])/2 + ifelse(runif(1,-1,1) > 0, 1, -1) * runif(1,0,d/3),
                 (coords[from, 2] + coords[to, 2])/2 + ifelse(runif(1,-1,1) > 0, 1, -1) * runif(1,0,d/3))
    
    # BEZIER LINES
    b = bezier(t=t <- seq(0, 1, length=100), p=rbind(coords[from, ], midpoint, coords[to, ]))
    
    # WE BUILD THE SPATIAL LINE OBJECT
    rowname = paste(from, to, sep="_")
    l = makeLineFromCoords(b, rowname)
    df = data.frame("from" = bati$C_PERCONST[from], "to" = bati$C_PERCONST[to], row.names=rowname) # WE GET THE PERIOD
    df$DIFF = abs(df$PER_CONST_FROM - df$PER_CONST_TO)
    l.df = SpatialLinesDataFrame(l, data=df)
    
    synapses[[i]] = l.df

  }
  
  res.df = do.call(rbind, synapses)
  writeOGR(res.df, "OUT", paste("mstCumul_perconst", periode), "ESRI Shapefile", overwrite=T) # ONE FILE PER PERIOD BUT CUMULATED
  
}