20 Nov

Get angle/orientation/direction of GPS track using R

Hey,

Today a bit of trSIGonometry with this function to compute the angle/orientation/direction between 2 points of a GPX file.

############################
#Function
############################
 
angle <- function(pt1, pt2) {
  conv = 180 / pi # conversion radians / degrees
  diff <- pt2 - pt1
  diff.abs <- abs(diff)
  if (diff[1]>=0 & diff[2]>=0){angle <- pi/2 - atan(diff.abs[2] / diff.abs[1])}
  if (diff[1]>=0 & diff[2]<=0){angle <- pi/2 + atan(diff.abs[2] / diff.abs[1])}
  if (diff[1]<=0 & diff[2]>=0){angle <- 3*pi/2 + atan(diff.abs[2] / diff.abs[1])}
  if (diff[1]<=0 & diff[2]<=0){angle <- pi + (pi/2- atan(diff.abs[2] / diff.abs[1]))}
  return(angle * conv)
}
 
############################
#Script
############################
 
library(rgdal)
library(raster)
library(data.table)
 
#Set folder
setwd("/home/GuruBlard/gpx/")
 
# Load file names
filenames <- list.files(pattern=".gpx$")
 
# Iteration through all gpx files
for (filename in filenames){
  # Load the GPX
 
  #ogrListLayers("export.gpx")
  gpx.trackpoints <- readOGR(dsn = filename,  layer= "track_points", stringsAsFactors = F)
 
  # Extract the coordinates as a matrix:
  gpx.tp <- coordinates(gpx.trackpoints)
 
  # #Loop through and add angle - the difference between the of the current point and the following one:
  for (tp in 1:(nrow(gpx.tp)-1)) {
    gpx.trackpoints$angle[tp] <- angle(gpx.tp[tp,], gpx.tp[tp+1,])
  }
 
  gpx.trackpoints@data$Date <- substr(as.character(gpx.trackpoints@data$time), 1,10)
  gpx.trackpoints@data$Hour <- substr(as.character(gpx.trackpoints@data$time),12,19)
  gpx.trackpoints@data$filename <- filename
  writeOGR(gpx.trackpoints, paste0(tools::file_path_sans_ext(filename), ".shp"), layer=tools::file_path_sans_ext(filename), driver="ESRI Shapefile")
}

Screenshot from 2017-11-20 14-53-09

13 Jul

Italic label in a plot using facet_wrap

Hi!
I just got asked by a journal to italicize some labels in a plot… So here is a small code snippet to assist you when you get asked the same:

factor1=rep(letters[1:3], each=3)
factor2=rep(1:3,times=3)
x=rep(1,9)
y=1:9
df=cbind.data.frame(factor1,factor2,x,y)
 
levels(df$factor1)= c("a"=expression(paste("factor_", italic("a"))),
                      "b"=expression(paste("factor_", italic("b"))),
                      "c"=expression(paste("factor_", italic("c"))))
 
ggplot(df, aes(x=x, y=x))+facet_grid(factor2~factor1, labeller=label_parsed)+geom_point()

And the result is:
italic

Thanks to Berengere Husson for the tip!

01 Jun

Bland-Altman diagram with ggplot2

Hello world!

Today, I’d like to share with you a code for plotting Bland-Altman diagrams.
The Bland-Altman diagram (Bland & Altman, 1986 and 1999), or difference plot, is a graphical method to compare two measurements techniques. In this graphical method the differences (or alternatively the ratios) between the two techniques are plotted against the averages of the two techniques. Alternatively (Krouwer, 2008) the differences can be plotted against one of the two methods, if this method is a reference or “gold standard” method.
Horizontal lines are drawn at the mean difference, and at the limits of agreement, which are defined as the mean difference plus and minus 1.96 times the standard deviation of the differences.

Bland.Altman.diagram <-
  function(x,y,id=FALSE, alpha = .05,rep.meas = FALSE,subject,idname = FALSE, xname = 'x',yname = 'y',...) {
 
    library(ggplot2)
	  library(RColorBrewer)
 
    #*** 1. Set a few constants
    z <- qnorm(1 - alpha / 2)  ## value of z corresponding to alpha
    d <- x - y               ## pair-wise differences
    m <- (x + y) / 2           ## pair-wise means
 
    #*** 2. Calculate mean difference
    d.mn <- mean(d,na.rm = TRUE)
 
    #*** 3. Calculate difference standard deviation
    if (rep.meas == FALSE) {
      d.sd = sqrt(var(d,na.rm = TRUE))
    } else {
      #*** 3a. Ensure subject is a factor variable
      if (!is.factor(subject))
        subject <- as.factor(subject)
 
      #*** 3b. Extract model information
      n <- length(levels(subject))      # Number of subjects
      model <- aov(d ~ subject)           # One way analysis of variance
      MSB <- anova(model)[[3]][1]       # Degrees of Freedom
      MSW <- anova(model)[[3]][2]       # Sums of Squares
 
      #*** 3c. Calculate number of complete pairs for each subject
      pairs <- NULL
      for (i in 1:length(levels(as.factor(subject)))) {
        pairs[i] <- sum(is.na(d[subject == levels(subject)[i]]) == FALSE)
      }
      Sig.dl <-
        (MSB - MSW) / ((sum(pairs) ^ 2 - sum(pairs ^ 2)) / ((n - 1) * sum(pairs)))
      d.sd <- sqrt(Sig.dl + MSW)
    }
 
    #*** 4. Calculate lower and upper confidence limits
    ucl <- d.mn + z * d.sd
    lcl <- d.mn - z * d.sd
    print(d.mn)
    print(ucl)
    print(lcl)
 
    #*** 5. Make Plot
    xlabstr = paste(c('Mean(',xname,', ',yname,')'), sep = "",collapse = "")
    ylabstr = paste(c(xname, ' - ', yname), sep = "",collapse = "")
    id = ifelse(id==FALSE, as.factor(rep(1,length(m))), as.factor(id))
    xy <- data.frame(m,d,id)
    xy$id <- as.factor(xy$id)
    maxm <- max(m)
 
 
    #scatterplot of x and y variables
    scatter <- ggplot(xy,aes(m, d)) 
 
 
    if(idname == FALSE){
      scatter <- scatter + geom_point(color='darkgrey') + scale_color_manual( guide=FALSE)
    } else {
      scatter <- scatter + geom_point(aes(color = id))+ scale_color_brewer(idname, palette="Set1") # Set1
    }
 
    scatter <- scatter + geom_abline(intercept = d.mn,slope = 0, colour = "black",size = 1,linetype = "dashed" ) +
      geom_abline(
        intercept = ucl,slope = 0,colour = "black",size = 1,linetype = "dotted"
      ) +
      geom_abline(
        intercept = lcl,slope = 0,colour = "black",size = 1,linetype = "dotted"
      ) +
      ylim(min(lcl, -max(xy$d), min(xy$d)), max(ucl, max(xy$d), -min(xy$d)))+
      xlab(xlabstr) +
      ylab(ylabstr) +
      #geom_text(x=maxm-2.4,y=d.mn,label='Mean',colour='black',size=4)+
      #geom_text(x=maxm-2.4,y=ucl,label='+1.96 SD',colour='black',size=4)+
      #geom_text(x=maxm-2.4,y=lcl,label='-1.96 SD',colour='black',size=4)+
      theme(
        plot.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(color="black", size = 1)
    )
 
    plot(scatter)
 
    values <- round(cbind(lcl,d.mn,ucl),4)
    colnames(values) <- c("LCL","Mean","UCL")
    if (rep.meas == FALSE){
      Output <- list(limits = values,Var = d.sd ^ 2)
    } else {
      Output <- list(limits = values,Var = Sig.dl,MSB = MSB,MSW = MSW)
    }
    return(Output)
  }
 
A <- c(-0.358, 0.788, 1.23, -0.338, -0.789, -0.255, 0.645, 0.506, 
       0.774, -0.511, -0.517, -0.391, 0.681, -2.037, 2.019, -0.447, 
       0.122, -0.412, 1.273, -2.165)
B <- c(0.121, 1.322, 1.929, -0.339, -0.515, -0.029, 1.322, 0.951, 
       0.799, -0.306, -0.158, 0.144, 1.132, -0.675, 2.534, -0.398, 0.537, 
       0.173, 1.508, -1.955)
 
Bland.Altman.diagram(A, B,alpha = .05,rep.meas = F,xname = 'A',yname = 'B')

The result is the following:
Rplot01

The ‘id’ field allows you to define different categories that will be colored differently on the graph:Rplot02

18 Apr

Online LaTeX diff tool

Hi PubliGISer,

You want to publish the revision of your awesome paper but you need a track change to meet reviewers’ requirements and you do not know how to get it? You are desperate and think that you will not be able to realize your dream of a successful research career?

Online LaTeX diff tool is the TOOL you are looking for.

Just copy paste the old and the new latex files and this tool will automatically compute the difference between the 2 documents.

Problem solved.

22 Nov

3D plot with a cut-off line

Did you ever notice that things IRL are in 3D? Well today, let’s do some 3D plots!

First, let’s define our three variables

library(rgl)
 
x <- 1:5/10
y <- 1:5
z <- x %o% y
z <- z + .2*z*runif(25) - .1*z

Now let’s suppose that we want to plot an additional line at a specific cut-off value of 0.6

cutoff <- 0.6

Let’s plot all that with the RGL package:

persp(x, y, z, theta=-35, phi=10,col="lightgrey",xlab="X factor", ylab="Why", zlab="Z-bra")
cLines <- contourLines(x,y,z,levels=c(cutoff))
lines(trans3d(x=cLines[[1]]$x, y=cLines[[1]]$y, z=cutoff,pmat=p ),col = 'red',lw=2,lt=2)

Here is the figure:
screenshot-from-2016-11-22-175501

22 Jan

Styling shapefiles for QGIS: creating symbology in R

Hello world!

Today, I’d like to share with y’all a function that I wrote to save symbology associated with a shapefile of points in R.
All you need is you shapefile in which you should have three columns:
1. col: indicating the color you want to associate with to your class in hex values (see the super http://www.color-hex.com/ website.).
2. label: the label of your class
3. value: the numeric value of your label

Enjoy!

Here is the function (you can also find it on my shiny GITHUB). You can also simply type in source_url(“https://rawgit.com/waldnerf/R/master/writeQML.R”) (! requires the devtools package).

writeQML <- function(inShp, col, label, value, filename){
  require(RColorBrewer)
  require(XML)
  inNames <- names(inShp)
  names(inShp)[which(names(inShp)==value)] <- "Value"
  names(inShp)[which(names(inShp)==col)] <- "v"
  names(inShp)[which(names(inShp)==label)] <- "Label"
 
  inShp$v <- unlist(lapply(as.character(inShp$v),function(x) paste(c(col2rgb(x),'255'),collapse = ',')))
  inShp$Label <- as.character(inShp$Label)
 
  inShp <- inShp[-which(duplicated(inShp@data[,c('Value','v','Label')])),]
 
  # set alpha
  alpha = 1
 
  base = newXMLNode("qgis")
  addAttributes(base,version="2.8.1-Wien", minimumScale="0", maximumScale="1e+08", simplifyDrawingHints="0", minLabelScale="0", maxLabelScale="1e+08",simplifyDrawingTol="1" ,simplifyMaxScale="1", hasScaleBasedVisibilityFlag="0", simplifyLocal="1", scaleBasedLabelVisibilityFlag="0")
  trans <- newXMLNode("transparencyLevelInt", 255)
  rend <- newXMLNode("renderer-v2", attrs = c(attr=value,symbollevels="0",type="categorizedSymbol"))
 
  # sort the categories
  categories <- newXMLNode("categories")
  category <- lapply(seq_along(inShp$Value),function(x){newXMLNode("category", attrs = c(symbol = as.character(x-1), value = inShp$Value[x], label = inShp$Label[x]))
  })
  addChildren(categories,category)
 
  # sort the symbols
  symbols <- newXMLNode("symbols")
  symbol <- lapply(seq_along(inShp$Value),function(x){dum.sym <- newXMLNode("symbol", attrs = c(outputUnit="MM",alpha=alpha,type="marker",name=as.character(x-1)))
  layer <- newXMLNode("layer", attrs =c(pass="0",class="SimpleMarker",locked="0"))
  prop <- newXMLNode("prop", attrs =c(k="color",v= inShp$v[x]))
  addChildren(layer, prop)
  addChildren(dum.sym, layer)
  }) 
  addChildren(symbols, symbol)
 
  # add categories and symbols to rend
  addChildren(rend, list(categories, symbols))
  addChildren(base, list(trans, rend))
 
  # save to qml-file
  writeLines(saveXML(base, prefix = "<!DOCTYPE qgis >"), filename)
}

And a self-sufficient example.

# Prepare the shapefile
inDf <- structure(list(x = c(25.582125, 25.880375, 24.955375, 26.066375, 25.409375, 25.902125),
                       y = c(-26.147625, -26.202625, -26.166125, -26.117625, -26.035875, -26.046375), 
                       Class = c(1L, 2L, 3L, 10L, 10L, 10L),
                       colT = c("#ffa500", "#aaeeff", "#a8c093", "#4584dc","#4584dc", "#4584dc"), 
                       attT = c("Stuff", "Other stuff", "Thing","Nice thing", "Nice thing", "Nice thing")), 
                       .Names = c("x", "y", "Class", "colT", "attT"),
                       row.names = c(NA, 6L), class = "data.frame")
 
spat.df <- SpatialPointsDataFrame(inDf[,c('x','y')], inDf)
projection(spat.df) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
shapefile(spat.df, filename='/to/your/path/example.shp',overwrite=T)
filename <- "/to/your/path/example.qml"
 
writeQML(spat.df, "colT", 'attT', 'Class', filename)
21 Jan

Dowload OpenStreetMap (OSM) data and read it in QGIS

Happy New GIS,

Be sure to get the last version of GDAL (Ubuntu installation)

Go to Geofabrik.
Download the compressed .osm.pbf files corresponding to our area (Continent, Country).

Using ogr to convert the pbf file into a SQlite database which can be used directly in QGIS with Add SpatiaLite Layer.

ogr2ogr -f "SQLite" -dsco SPATIALITE=YES -spat xmin ymin xmax ymax GuruBlard_Region.db OSM_continent.pbf

where xmin ymin xmax ymax:
spatial query extents, in the SRS of the source layer(s) (or the one specified with -spat_srs). Only features whose geometry intersects the extents will be selected. The geometries will not be clipped unless -clipsrc is specified

For further information, visit this post of Anita Graser.

08 Jun

Assemble png into a grid using ImageMagick and R

Holà GISpa,

Here is a R function that uses the montage function of ImageBrick to assemble png in one grid. All the parameters are explained here.

pattern='CoolestPngEver_'
In<-list.files("/export/homes/gurublard/", full.names=T, pattern=pattern)
Out<-paste0(/export/homes/gurublard/", pattern,"all.png")
command<-paste0("montage -label '%f' -density 300 -tile 4x0 -geometry +5+50 -border 10 ", paste(In, collapse=" "), " ", Out)
system(command)

ImageMagick

27 Mar

Calculate hillshade with gdaldem based on the elevation and azimuth angle calculated from the acquisition time

Hi GIS fellows ! Here is a spring present, let’s see how to simply create an hillshade based on the acquisition date and time !

sunPosition <- function(year, month, day, hour=12, min=0, sec=0,
                        lat=46.5, long=6.5) {
 
  twopi <- 2 * pi
  deg2rad <- pi / 180
 
  # Get day of the year, e.g. Feb 1 = 32, Mar 1 = 61 on leap years
  month.days <- c(0,31,28,31,30,31,30,31,31,30,31,30)
  day <- day + cumsum(month.days)[month]
  leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & 
    day >= 60 & !(month==2 & day==60)
  day[leapdays] <- day[leapdays] + 1
 
  # Get Julian date - 2400000
  hour <- hour + min / 60 + sec / 3600 # hour plus fraction
  delta <- year - 1949
  leap <- trunc(delta / 4) # former leapyears
  jd <- 32916.5 + delta * 365 + leap + day + hour / 24
 
  # The input to the Atronomer's almanach is the difference between
  # the Julian date and JD 2451545.0 (noon, 1 January 2000)
  time <- jd - 51545.
 
  # Ecliptic coordinates
 
  # Mean longitude
  mnlong <- 280.460 + .9856474 * time
  mnlong <- mnlong %% 360
  mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360
 
  # Mean anomaly
  mnanom <- 357.528 + .9856003 * time
  mnanom <- mnanom %% 360
  mnanom[mnanom < 0] <- mnanom[mnanom < 0] + 360
  mnanom <- mnanom * deg2rad
 
  # Ecliptic longitude and obliquity of ecliptic
  eclong <- mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)
  eclong <- eclong %% 360
  eclong[eclong < 0] <- eclong[eclong < 0] + 360
  oblqec <- 23.439 - 0.0000004 * time
  eclong <- eclong * deg2rad
  oblqec <- oblqec * deg2rad
 
  # Celestial coordinates
  # Right ascension and declination
  num <- cos(oblqec) * sin(eclong)
  den <- cos(eclong)
  ra <- atan(num / den)
  ra[den < 0] <- ra[den < 0] + pi
  ra[den >= 0 & num < 0] <- ra[den >= 0 & num < 0] + twopi
  dec <- asin(sin(oblqec) * sin(eclong))
 
  # Local coordinates
  # Greenwich mean sidereal time
  gmst <- 6.697375 + .0657098242 * time + hour
  gmst <- gmst %% 24
  gmst[gmst < 0] <- gmst[gmst < 0] + 24.
 
  # Local mean sidereal time
  lmst <- gmst + long / 15.
  lmst <- lmst %% 24.
  lmst[lmst < 0] <- lmst[lmst < 0] + 24.
  lmst <- lmst * 15. * deg2rad
 
  # Hour angle
  ha <- lmst - ra
  ha[ha < -pi] <- ha[ha < -pi] + twopi
  ha[ha > pi] <- ha[ha > pi] - twopi
 
  # Latitude to radians
  lat <- lat * deg2rad
 
  # Azimuth and elevation
  el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
  az <- asin(-cos(dec) * sin(ha) / cos(el))
 
  # For logic and names, see Spencer, J.W. 1989. Solar Energy. 42(4):353
  cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
  sinAzNeg <- (sin(az) < 0)
  az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
  az[!cosAzPos] <- pi - az[!cosAzPos]
 
 
  el <- el / deg2rad
  az <- az / deg2rad
  lat <- lat / deg2rad
 
  return(list(elevation=el, azimuth=az))
}
 
 
 
# calculate the sun position from date and time
position<-sunPosition(2012, 05, 25, hour=15, min=40, sec=42.460,lat=d_wgs@coords[2], long=d_wgs@coords[1])
 
system(paste('gdaldem hillshade',
             '-of GTiff',
             '-az',position$azimuth,
             '-alt',position$elevation,
             paste('/path_to_data/dtm_file','.tif',sep=''),
             paste('/path_to_data/dtm_file','_hillshade.tif',sep='')))