View source for Main Page/Research/MSB/processing model/20090607/line edit to points.R

<pre>
# line to point from PostGIS
source("~/tools/dump.postgis.z.coords.R")
source("~/tools/calc.distance.R")
source("http://gis.washington.edu/phurvitz/R/functions.R")

line_edit_to_points <- function(sid) {
    
    # start time 
    time.start <- Sys.time()
    cat(paste("Started at", time.start, "\n"))
    # subject dir
    if (.Platform$OS.type=="windows") {
        basedir <- "q:/msb/"
    } else {
        basedir <- "/home/phurvitz/msb/"
    }
    subject.dir <- paste(basedir, sid, sep="")   
    setwd(subject.dir) 
    
    # convert the *_line_edit lines to points (dump.postgis.z.coords.R)
    input.table.name <- sprintf("%s_line_edit", sid)
    seq_id.trip_id <- dump.postgis.z.coords(sid, input.table.name)
    
    # read in original point data
    cat("reading in original data ....\n")
    input.point.file <- sprintf("%s_points.csv", sid)
    dat.raw <- read.csv(input.point.file, as.is=T)
    
    # join the original data onto the edited point data using the seq_id
    cat("merging....\n")
    dat.joined <- merge(seq_id.trip_id, dat.raw, by="seq_id")
    dat.joined$bout <- NULL
    colnames(dat.joined) <- fix.colnames(dat.joined)
    
    # distance
    dat.joined <- with(dat.joined, calc.distance(dat.joined, x.waspn, y.waspn, strpad="   "))
    
    # group by trip_id
    trips <- split(dat.joined, dat.joined$trip.id)
    trips.joined <- NULL
    
    # for each trip
    for (trip.num in 1:length(trips)) {
        trip <- trips[[trip.num]]
        movetype <- trip$movetype[1]
        cat(paste(trip.num, "of", length(trips), "calculating straightened points....\n"))
        # cumulative distance
        cumdist <- cumsum(trip$dist.next)
        cumdist <- c(0, cumdist[-length(cumdist)])
        # length 
        length.ft <- tail(cumdist,1)
        # duration
        duration.s <- with(trip, msg41.navutc.time.t[nrow(trip)] - msg41.navutc.time.t[1])
        # velocity
        trip$vel.fps <- length.ft/duration.s
        x0 <- trip[1,"x.waspn"]
        y0 <- trip[1,"y.waspn"]
        trip$x.straight <- x0 + cumdist
        trip$y.straight <- y0
        

        # radius counts
        x.straight <- trip$x.straight
        for (radius in c(50, 75, 100, 125, 150)) {
            cat(paste("counting points in radius", radius, "....\n"))
            # moving only
            if (movetype=="m") {
                pointcount <- NULL
                for (i in 1:length(x.straight)) {
                    x.index <- x.straight[i]
                    x.start <- x.index - radius
                    x.end <- x.index + radius
                    pointcount <- c(pointcount, length(x.straight[x.straight >= x.start & x.straight < x.end]))
                }
                
                cmd <- sprintf("trip$rcount.%s <- pointcount", radius)
                eval(parse(text=cmd))
            } else {
                cmd <- sprintf("trip$rcount.%s <- -9", radius)
                eval(parse(text=cmd))            
            }
        }

        trips.joined <- rbind(trips.joined, trip)
    }
    
    # write out a csv
    scipen <- options()$scipen
    options(scipen=10)
    colnames(trips.joined) <- unfix.colnames(trips.joined)
    out.csv <- sprintf("%s_trips.joined.csv", sid)
    write.csv(trips.joined, file=out.csv, row.names=F)
    options(scipen=scipen)
    
    # write the SQL to this into PostGIS
    outsql <- "trips.joined.sql"
    out.table.name <- sprintf("%s_trip_points", sid)
    cat(file=outsql, 
    sprintf("DROP TABLE IF EXISTS %s CASCADE;", out.table.name),
    sprintf("CREATE TABLE %s (", out.table.name),
    "seq_id integer,", 
    "trip_id integer,", 
    "movetype varchar(1),", 
    "sid varchar(3),", 
    "fileid integer,", 
    "id integer,", 
    "packetindex integer,", 
    "repaired_msg41_navutc_time_t integer,", 
    "elevationdata varchar(1),", 
    "movementdata varchar(1),", 
    "msg41_longitude double precision,", 
    "msg41_latitude double precision,", 
    "msg41_altitude_from_ellipsoid double precision,", 
    "msg41_navutc_time_t integer,", 
    "msb_timestamp integer,", 
    "repaired_msb_timestamp integer,", 
    "uwar_timestamp_offset integer,", 
    "repaired_uwar_timestamp_offset integer,", 
    "msg41_est_horz_pos_err real,", 
    "msg41_est_vert_pos_err real,", 
    "msg41_climb_rate real,", 
    "msg41_number_of_svs_in_fix integer,", 
    "msg41_hdop real,", 
    "gps_utctim timestamp,", 
    "gps_loctim timestamp,", 
    "allmove varchar(1),", 
    "allmove_num integer,", 
    "x_waspn double precision,", 
    "y_waspn double precision,", 
    "dist_prev real,", 
    "dist_next real,", 
    "dist_prevnext real,", 
    "vel_fps real,",
    "x_straight double precision,", 
    "y_straight double precision,",
    sep="\n")
    
    # a field for each radius
    radii <- c(50, 75, 100, 125, 150)
    for (i in 1:length(radii) ) {
        radius <- radii[i]
        if (i < length(radii)) {        
            field <- sprintf("rcount_%s", radius)
            cat(file=outsql, append=T, sprintf("%s integer,\n", field))
        } else {
            field <- sprintf("rcount_%s", radius)
            cat(file=outsql, append=T, sprintf("%s integer);\n", field))        
        }
    }
    
    cat(file=outsql, append=T,
    sprintf("COPY %s FROM '%s/%s' CSV HEADER;", out.table.name, subject.dir, out.csv),
    sprintf("SELECT AddGEometryColumn('', '%s', 'point_waspn', 2926, 'POINT', 2);", out.table.name),
    sprintf("SELECT AddGEometryColumn('', '%s', 'point_waspn_straight', 2926, 'POINT', 2);", out.table.name),
    sprintf("UPDATE %s SET point_waspn = PointFromText('POINT(' || x_waspn || ' ' || y_waspn || ')', 2926);", out.table.name),
    sprintf("UPDATE %s SET point_waspn_straight = PointFromText('POINT(' || x_straight || ' ' || y_straight || ')', 2926);", out.table.name),
    sep="\n")
    
    # run the SQL
    system(sprintf("psql msb -f %s", outsql))
    
    # column names 
    colnames(trips.joined) <- fix.colnames(trips.joined)
    
    # dump points
    cmd.dump <- sprintf("pgsql2shp -f %s_edited_points.shp -g point_waspn msb %s", sid, out.table.name)
    system(cmd.dump)
    cmd.dump <- sprintf("pgsql2shp -f %s_edited_points_straight.shp -g point_waspn_straight msb \"SELECT * FROM %s WHERE movetype = 'm'\"", sid, out.table.name)
    system(cmd.dump)
    sprintf("Created %s_edited_points.shp and %s_edited_points_straight.shp\n", sid, sid)

    # start time 
    time.end <- Sys.time()
    # duration
    duration <- round(difftime(time.end, time.start, units="s"), 0)
    cat(paste("Ended at", time.end, "\n"))

    cat(paste("Processing time:", duration, "s\n"))
    
    return(trips.joined)
    
}
</pre>