Skip to content

Commit

Permalink
Various improvements to fit code, more comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@lordsutch
lordsutch committed Dec 11, 2012
1 parent c6ff332 commit 28246c4
Show file tree
Hide file tree
Showing 2 changed files with 175 additions and 51 deletions.
153 changes: 120 additions & 33 deletions makeroads.R
View file
Expand Up @@ -15,14 +15,15 @@
## terms of the GNU General Public License, version 2 or later.

library(rgeos)
library(princurve)
library(rgdal)
library(maptools)
library(geosphere)
##library(LPCM)
library(princurve)

parseGPXfile <- function(filename) {
timesep <- 60 ## Split if points are more than 'timesep' seconds apart
parseGPXfile <- function(filename, tracksfrompoints=FALSE) {
timesep <- 60 ## Split if points are more than 'timesep' seconds apart
distsep <- 100 ## Only make tracks if points less than this distance apart

converted <- try(readOGR(filename, 'track_points'), TRUE)

Expand All @@ -34,6 +35,31 @@ parseGPXfile <- function(filename) {
points <- coordinates(converted)[!trackrows,,drop=FALSE]
lines <- NULL

## Try to assemble anonymized points into lines
if(tracksfrompoints) {
start <- 1
linelist <- NULL
npoints <- NULL

for(j in seq_int(nrow(points)-1)) {
dist <- distHaversine(points[j,], points[j+1,])
if(dist > distsep) {
seg <- points[start:j,,drop=FALSE]
if(nrow(seg) > 1) {
linelist <- c(linelist, Line(seg))
} else {
npoints <- rbind(npoints, seg)
}
start <- j+1
}
}

if(!is.null(linelist)) {
lines <- c(lines, Lines(linelist, i))
}
points <- npoints
}

i <- 1
alltracks <- converted[trackrows,]
for(tnum in unique(alltracks$track_fid)) {
Expand All @@ -51,7 +77,7 @@ parseGPXfile <- function(filename) {
linelist <- NULL
for(seg in splits) {
segpoints <- coordinates(thisseg)[start:seg,,drop=FALSE]
str(segpoints)
##str(segpoints)
if(!is.null(segpoints)) {
if(nrow(segpoints) > 1) {
linelist <- c(linelist, Line(segpoints))
Expand Down Expand Up @@ -79,6 +105,9 @@ parseGPXfile <- function(filename) {

tracks <- NULL
spoints <- NULL

str(points)

if(!is.null(lines))
tracks <- SpatialLines(lines, proj4string=CRS(proj))
if(!is.null(points))
Expand All @@ -97,7 +126,7 @@ fetchGPXpage <- function(left, bottom, right, top, page=0) {
print(filename1)
download.file(pageurl, filename1, cacheOK=FALSE)

converted <- parseGPXfile(filename1)
converted <- parseGPXfile(filename1, tracksfrompoints=closetracks)

if(!debug) unlink(filename1)
converted
Expand Down Expand Up @@ -149,7 +178,7 @@ getOSMtracksFiles <- function(...) {
tracks <- NULL
points <- NULL
for(filename in filenames) {
converted <- parseGPXfile(filename)
converted <- parseGPXfile(filename, tracksfrompoints=closetracks)
if(!is.null(converted)) {
if(is.null(tracks)) {
tracks <- converted$tracks
Expand Down Expand Up @@ -193,28 +222,40 @@ splitTracksAngle <- function(tracks, splitangle, splitdistance) {
az <- trackAzimuth(cds, type='abdali')
if(length(az) > 3) {
laz <- c(tail(az, -3), NA, NA, NA)
show(length(az))
show(length(laz))
ldiff <- pmin((az-laz) %% 360, (laz-az) %% 360)
ldiff[is.na(ldiff)] <- 0
str(ldiff)
} else {
ldiff <- rep(180, length(az))
ldiff <- rep(0, length(az))
}
repeat {
##str(d)
firstangle <- az[1]

adiff <- pmin((az-firstangle) %% 360, (firstangle-az) %% 360, ldiff)
distfromstart <- cumsum(d)
str(distfromstart)
pos <- match(TRUE, adiff > splitangle &
distfromstart > minanglesplit)
if(is.na(pos)) {
linelist <- c(linelist, Line(pcds))
adiff <- c(0,pmin((az-firstangle) %% 360, (firstangle-az) %% 360))
adiff <- pmax(adiff, c(0,ldiff))

distfromstart <- c(0, cumsum(d))
##str(distfromstart)
x <- which(adiff > splitangle & distfromstart > minanglesplit)
str(x)
if(length(x) > 0 && x[[1]] == 1) {
x <- x[-1]
}
if(length(x) < 1) {
if(nrow(pcds) > 1)
linelist <- c(linelist, Line(pcds))
else
points <- rbind(points, pcds)
break
}
str(pos)
linelist <- c(linelist, Line(pcds[1:pos,,drop=FALSE]))
pos <- x[1]
bit <- pcds[1:pos,,drop=FALSE]
if(nrow(bit) > 1) {
linelist <- c(linelist, Line(bit))
} else {
points <- rbind(points, bit)
}
cds <- cds[-(1:pos-1),,drop=FALSE]
pcds <- pcds[-(1:pos-1),,drop=FALSE]
az <- az[-(1:pos-1),drop=FALSE]
Expand Down Expand Up @@ -268,7 +309,7 @@ findClosestTrackToPoint <- function(point, tracks, tracklist) {
closestTrack <- which.min(dists)
show(closestTrack)

if(mindist > maxtrackdist)
if(mindist > maxpointdist)
return(0)

for(t in seq_along(tracklist)) {
Expand Down Expand Up @@ -310,16 +351,16 @@ consolidateTracks <- function(tracks, points) {
ylim <- c(bounds['y','min'], bounds['y','max'])
xlim <- c(bounds['x','min'], bounds['x','max'])

newtracks <- flattenSpatialLines(tracks)
##newtracks <- flattenSpatialLines(tracks)

## Sort by length
newtracks <- sortTracks(newtracks)
##newtracks <- sortTracks(newtracks)
newtracks <- tracks

bearings <- integer(length(newtracks))
for(i in seq_along(newtracks)) {
bearings[[i]] <- calcbearings(newtracks[i])
}
loosepoints <- NULL

t <- 2
while(t <= length(newtracks)) {
Expand Down Expand Up @@ -378,9 +419,8 @@ consolidateTracks <- function(tracks, points) {
closeenough <- FALSE

for(ctrack in tracklist[[v]]) {
strack <- gSimplify(newtracks[ctrack], 1, TRUE)
buffzone <- gBuffer(strack, width=maxtrackdist)
fatzone <- gBuffer(strack, width=separation)
buffzone <- gBuffer(newtracks[ctrack], width=maxtrackdist)
fatzone <- gBuffer(newtracks[ctrack], width=separation)

if(!gCrosses(track, buffzone))
next
Expand All @@ -393,8 +433,8 @@ consolidateTracks <- function(tracks, points) {
insidepart <- gIntersection(track, fatzone)
outsidepart <- gDifference(track, fatzone)

str(insidepart)
str(outsidepart)
##str(insidepart)
##str(outsidepart)

lines(newtracks[ctrack], col='green')
if(!is.null(outsidepart))
Expand Down Expand Up @@ -537,17 +577,25 @@ showtracks <- function(tracks, basetrack, others) {
}

sdistances <- function(curve, track, projection) {
dvec <- numeric(nrow(track))
ptrack <- project(as.matrix(track), projection, inv=TRUE)
ctrack <- project(as.matrix(curve$s), projection, inv=TRUE)
ctrack <- project(as.matrix(curve), projection, inv=TRUE)

for(r in 1:nrow(track)) {
dvec[[r]] <- distHaversine(ptrack[r,], ctrack[r,])
}
##str(ptrack)
##str(ctrack)

dvec <- distHaversine(ptrack, ctrack)
str(dvec)
str(sd(dvec))

sdists <- dvec/sd(dvec)
sdists
}

setupForFit <- function(points, tracklist, trackforpoints, tcount) {
thesepoints <- points[trackforpoints == tcount]
coordinates(thesepoints)
}

fitpcurve <- function(track, projection) {
olen <- nrow(track)
bandwidth <- round(min(0.1, 30/olen), 3)
Expand All @@ -573,7 +621,46 @@ fitpcurve <- function(track, projection) {
} else {
show('Empty refit?')
}
curve
curve[curve$tag,]
}

## Alternative using local principal curves algorithm. At present,
## more fiddly than the loess fit, so disabled.

## Could weight using hdop if available...
fitpcurve.lpc <- function(track, projection, bandwidth=0.14) {
track <- as.matrix(track)
olen <- nrow(track)
show(paste('Fitting curve with', olen, 'points; bandwidth', bandwidth))

x0 <- 0 #track[sample(olen, 10),]
#str(x0)
spoints <- max(floor(olen/50), 5)

curve <- lpc(track, h=bandwidth, scaled=TRUE, pen=3, depth=3, x0=x0,
control=lpc.control(mult=spoints))
##plot(curve)
ucurve <- lpc.spline(curve, project=TRUE)
ucurve <- unscale(ucurve)
str(ucurve)

d <- sdistances(ucurve$closest.coords, track, projection)
show(which(d >= 3))
weights <- pmin(pmax((4-d), 0.2), 1)
str(d)
str(weights)

if(FALSE) { # any(weights) < 1
show(paste('Refitting weighted curve'))

curve <- lpc(track, h=bandwidth, scaled=TRUE, weights=weights,
depth=3, pen=3, x0=x0,
control=lpc.control(mult=spoints))
##plot(curve)
ucurve <- lpc.spline(curve, project=TRUE)
ucurve <- unscale(ucurve)
}
ucurve
}

gpsbabel.out <- function(infile, outfile) {
Expand Down

0 comments on commit 28246c4

Please sign in to comment.