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

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:

10 Nov

Count points in polygons

Hi guys,

Here is a quick r-snippet to count points inside polygons:

x <- getData('GADM', country='ITA', level=1)
# [1] "SpatialPolygonsDataFrame"
# attr(,"package")
# [1] "sp"
# sample random points
p <- spsample(x, n=300, type="random")
p <- SpatialPointsDataFrame(p, data.frame(id=1:300))
# [1] " +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
# [1] " +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
plot(p, col="red" , add=TRUE)

Here is the figure:

res <- over(p, x)
table(res$NAME_1) # count points
#               Abruzzo                Apulia            Basilicata
#                    11                    20                     9
#              Calabria              Campania        Emilia-Romagna
#                    16                     8                    25
# Friuli-Venezia Giulia                 Lazio               Liguria
#                     7                    14                     7
#             Lombardia                Marche                Molise
#                    22                     4                     3
#              Piemonte              Sardegna                Sicily
#                    35                    18                    21
#               Toscana   Trentino-Alto Adige                Umbria
#                    33                    15                     6
#         Valle d'Aosta                Veneto
#                     4                    22
15 May

Load googlesheet data in R

Dear R-addict,

At the era of cloud-based computing, functions such as read.csv, read.xls or even fread are really old-fashioned.
Here is a chunk of code that will allow you to load data from a googlesheet:

# Download a sheet
# Download a sheet
url <- 'https://docs.google.com/spreadsheets/d/1XBs7p44-djCPmN4TnPEgboUVAdB2mChbAlCjqnVOyQ0'
a <- gsheet2tbl(url)
12 May

Download daymet data with R

Hey folks ! It’s been a while !

Interested in meteorological data over US ? Here is a way to easily download Daymet data inspired from DaymetR package.

Daymet is a collection of algorithms and computer software designed to interpolate and extrapolate from daily meteorological observations to produce gridded estimates of daily weather parameters. Weather parameters generated include daily surfaces of minimum and maximum temperature, precipitation, humidity, and radiation produced on a 1 km x 1 km gridded surface over the conterminous United States, Mexico, and Southern Canada.

Download the “tiles” shapefile here.

param = c("vp", "tmin", "tmax", "swe", "srad", "prcp", "dayl") %see here https://daymet.ornl.gov/overview.html for the available variable
tiles=c(11922:11925, 11742:11745) %id of the tiles of interest (cfr shapefile)
for (i in year_range) {
  for (j in tiles) {
    for (k in param) {
      download_string = sprintf("https://daymet.ornl.gov/thredds/fileServer/ornldaac/1219/tiles/%s/%s_%s/%s.nc", 
                                i, j, i, k)
      daymet_file = paste(k, "_", i, "_", j, ".nc", 
                          sep = "")
      cat(paste("Downloading DAYMET data for tile: ", 
                j, "; year: ", i, "; product: ", k, "\n", 
                sep = ""))
      try(download(download_string, daymet_file, quiet = TRUE, 
                   mode = "wb"), silent = FALSE)
10 Mar

Efficient manuscript revision with latexdiff

Hi there,

Here is a nice tip to highlight the revisions in your manuscripts before resubmitting to your favorite journals: use latexdiff!

latexdiff --type=CTRADITIONAL old_shitty_manuscript.tex new_shiny_manuscript.tex > diff.tex

You can then easily customize how differences appear in your pdf by changing the commands in the preamble.
Here are my settings:

\providecommand{\DIFadd}[1]{{\protect\color{blue} #1}} %DIF PREAMBLE
\providecommand{\DIFdel}[1]{} %DIF PREAMBLE

New text appears in blue and old text is simply removed.

19 Feb

Color a raster with gdal using gdaldem

Hi Gigis,
Have you never look for a simple way to color a raster without manualy modifying the symbology or use a python/r script to build a RAT?
Here is the solution for lazy gis people like me 😉
My example is a raster Geotif with 0 for land, 1 for water and 255 for no data. Most of the time, when you open it in a GIS, you have black overview because the contrast is adjusted on the 255 value. Here is how to proceed:

1/ Write a simple text file colortable.txt such as:
0 white
1 blue
255 grey
2/gdaldem color-relief -of VRT sourcefile.tif colortable.txt out.vrt destfile.tif

3/ Just open the VRT, and the magic is there!

Khan-Guru, the Raster Magician

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


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){
  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]))
  # 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.

31 Aug

Merging two dataframe while keeping the original rows’ order

Hello GIStouille,

Here is a nice function to merge two dataframe while keeping the original rows’ order (source).
keep_order=1 keeps order of the first dataframe (x).
keep_order=2 keeps order of the second dataframe (y).

merge.with.order <- function(x,y, ..., sort = T, keep_order)
  # this function works just like merge, only that it adds the option to return the merged data.frame ordered by x (1) or by y (2)
  add.id.column.to.data <- function(DATA)
    data.frame(DATA, id... = seq_len(nrow(DATA)))
  # add.id.column.to.data(data.frame(x = rnorm(5), x2 = rnorm(5)))
  order.by.id...and.remove.it <- function(DATA)
    # gets in a data.frame with the "id..." column.  Orders by it and returns it
    if(!any(colnames(DATA)=="id...")) stop("The function order.by.id...and.remove.it only works with data.frame objects which includes the 'id...' order column")
    ss_r <- order(DATA$id...)
    ss_c <- colnames(DATA) != "id..."
    DATA[ss_r, ss_c]
  # tmp <- function(x) x==1; 1	# why we must check what to do if it is missing or not...
  # tmp()
    if(keep_order == 1) return(order.by.id...and.remove.it(merge(x=add.id.column.to.data(x),y=y,..., sort = FALSE)))
    if(keep_order == 2) return(order.by.id...and.remove.it(merge(x=x,y=add.id.column.to.data(y),..., sort = FALSE)))
    # if you didn't get "return" by now - issue a warning.
    warning("The function merge.with.order only accepts NULL/1/2 values for the keep_order variable")
  } else {return(merge(x=x,y=y,..., sort = sort))}
25 Aug

How to convert hdf files in tif?

Hi crazyGISfans,

Hdf format may be sometimes difficult to manage. We prefer working with geotif.
How do you do convert hdf in tif in R with gdal_translate?

files<-list.files(, pattern=".hdf$")
for (file in files){
  gdal_translate(src_dataset=get_subdatasets(file)[band], dst_dataset=paste0(file_path_sans_ext(file),"_B", band,'.tif'))

That’s it.