code for multilateration of cell towers in R

plotcells.R

library(ggmap)
library(plyr)
library(sp)
library(geosphere)

library(doMC)
registerDoMC()

allcells <- data.frame()
for(cell in list.files(pattern="^cellinfolte.*[.]csv$")) {
    newcells <- read.csv(cell)
    allcells <- rbind.fill(allcells, newcells)
}

allcells <- subset(allcells, select=-c(timestamp))
allcells <- subset(allcells, !is.na(allcells$timingAdvance) & !is.na(allcells$gci) & allcells$gci != "")

allcells$sector <- strtoi(substr(allcells$gci, 7, 8), 16)
allcells$base <- factor(substr(allcells$gci, 2, 6), ordered=T)

## Band 25 and 26 on same towers - this only works unmodified in Sprint GA market
allcells$baseDec <- strtoi(allcells$base, 16L)
allcells$base2 <- substr(allcells$gci, 2, 6)

## allcells$base2[allcells$band == 26] <- toupper(as.character(as.hexmode(477640-allcells$baseDec[allcells$band == 26])))
allcells$base2[allcells$baseDec >= 0x3A500] <- toupper(as.character(as.hexmode(477640-allcells$baseDec[allcells$baseDec >= 0x3A500])))

allcells$base2 <- factor(allcells$base2, ordered=T)

allcells <- unique(allcells)

allcells$distance <- allcells$timingAdvance*159.85
allcells$estDistance <- pmax(0, -2620.111-52.208*allcells$rsrp)

coords = cbind(allcells$longitude, allcells$latitude)
cellpoints <- SpatialPoints(coords, proj4string=CRS("+init=epsg:4326"))
area <- bbox(cellpoints)

basemap <- get_map(location=area, source="google", zoom=11)

## http://gis.stackexchange.com/questions/93126/trilateration-algorithm-for-n-amount-of-points-in-r?noredirect=1&lq=1
multilaterate <- function(x) {
    meanpos <- geomean(x[c('longitude', 'latitude')])
    if(length(x$latitude) < 3) {
        return(data.frame(latitude=meanpos[2], longitude=meanpos[1],
                          points=length(x$latitude)))
    }

    ch <- which.min(x$distance)
    sec <- x$sector[ch]

    a <- 0
    b <- 360
    ## Try to guess based on sector
    if(sec %in% c(1, 25)) {
        a <- 120
        b <- 240
    } else if (sec %in% c(2, 26)) {
        a <- 240
        b <- 360
    } else if (sec %in% c(3, 27)) {
        a <- 0
        b <- 120
    }
    
    startpoint <- destPoint(x[ch, c('longitude', 'latitude')],
                            runif(1, a, b),
                            x$distance[ch])

    n <- nls(distance ~ distm(data.frame(longitude, latitude),
                              c(lng_solution, lat_solution)),
             data=x, trace=TRUE, control=nls.control(warnOnly=TRUE),
             start=list(lng_solution=startpoint[1],
                        lat_solution=startpoint[2]))
    ##print(summary(n))
    data.frame(latitude=coef(n)[2], longitude=coef(n)[1],
               points=length(x$latitude))
}
centroids <- ddply(allcells, .(base2), multilaterate, .parallel=TRUE)
print(centroids)

drawnmap2 <- ggmap(basemap)
drawnmap2 <- drawnmap2 + geom_point(data=allcells, aes(x=longitude, y=latitude, colour=base2), alpha=.05, size=2) + theme(legend.position="none")
drawnmap2 <- drawnmap2 + geom_text(data=centroids, aes(x=longitude, y=latitude, label=base2, colour=base2), size=5)
drawnmap2

drawnmap1 <- ggmap(basemap)
drawnmap1 <- drawnmap1 + geom_point(data=allcells, aes(x=longitude, y=latitude, colour=as.factor(sector)), alpha=.05, size=2) + guides(colour = guide_legend(override.aes = list(alpha = 1)))
drawnmap1 <- drawnmap1 + geom_text(data=centroids, aes(x=longitude, y=latitude, label=base2), size=5)
drawnmap1

centroids <- ddply(allcells, .(base), multilaterate, .parallel=TRUE)
print(centroids)

drawnmap3 <- ggmap(basemap)
drawnmap3 <- drawnmap3 + geom_point(data=allcells, aes(x=longitude, y=latitude, colour=base), alpha=.05, size=2) + theme(legend.position="none")
drawnmap3 <- drawnmap3 + geom_text(data=centroids, aes(x=longitude, y=latitude, label=base, colour=base), size=5)
drawnmap3

pdf(file="Rplots.pdf")
drawnmap1
drawnmap2
drawnmap3
dev.off()