19 Feb

Larger and smaller extent shared by different rasters

Hello GISton,

Here is a way to get the larger (red) and the smaller (green) extent shared by different rasters.


extent_raster

library(raster)
 
# dummy extent from your rasters, instead use lapply(raster list, extent)
a <- raster(nrows=884, ncols=804, xmn=-45.85728, xmx=-43.76855, ymn=-2.388705, ymx=-0.5181549)
b <- raster(nrows=884, ncols=804, xmn=-45.87077, xmx=-43.78204, ymn=-2.388727, ymx=-0.5208711)
c <- raster(nrows=884, ncols=804, xmn=-45.81952, xmx=-43.7173,  ymn=-2.405129, ymx=-0.5154312)
 
a[] <- 1
b[] <- 2
c[] <- 3
 
plot(a, xlim=c(-45.8,-43.7), ylim=c(-2.41, -0.5))
par(new=TRUE)
plot(b, xlim=c(-45.8,-43.7), ylim=c(-2.41, -0.5))
 
a<-extent(-45, -30, -20, -10)
b<-extent(-55, -35, -25, -5) 
c<-extent(-40 ,-25 , -15 ,0)
extent_list<-list(a, b, c)
 
# make a matrix out of it, each column represents a raster, rows the values
extent_list<-lapply(extent_list, as.matrix)
matrix_extent<-matrix(unlist(extent_list), ncol=length(extent_list))
rownames(matrix_extent)<-c("xmin", "ymin", "xmax", "ymax")
 
# create an extent that covers all the individual extents
larger_extent<-extent(min(matrix_extent[1,]), max(matrix_extent[3,]),
                      min(matrix_extent[2,]), max(matrix_extent[4,]))
 
# create the larger extent shared by all the individual extents
smaller_extent<-extent(max(matrix_extent[1,]), min(matrix_extent[3,]),
                       max(matrix_extent[2,]), min(matrix_extent[4,]))
 
# Plot results
 
a.shp<-as(a, 'SpatialPolygons')
b.shp<-as(b, 'SpatialPolygons')
c.shp<-as(c, 'SpatialPolygons')
larger_extent.shp<-as(larger_extent, 'SpatialPolygons')
smaller_extent.shp<-as(smaller_extent, 'SpatialPolygons')
 
plot(a.shp,  xlim=c(larger_extent@xmin, larger_extent@xmax), ylim=c(larger_extent@ymin, larger_extent@ymax))
plot(b.shp,  xlim=c(larger_extent@xmin, larger_extent@xmax), ylim=c(larger_extent@ymin, larger_extent@ymax), add=TRUE)
plot(c.shp,  xlim=c(larger_extent@xmin, larger_extent@xmax), ylim=c(larger_extent@ymin, larger_extent@ymax), add=TRUE)
plot(larger_extent.shp,  xlim=c(larger_extent@xmin, larger_extent@xmax), ylim=c(larger_extent@ymin, larger_extent@ymax), 
     add=TRUE, lwd=3, border="red")
plot(smaller_extent.shp,  xlim=c(larger_extent@xmin, larger_extent@xmax), ylim=c(larger_extent@ymin, larger_extent@ymax), 
     add=TRUE, lwd=3, border="green")
13 Jan

Sample raster easily

Hey GISette! What’s up?

Here is a simple way to sample a raster as Raster object or SpatialPointsDataFrame using the library raster in r.

3plot

library(raster)
r <- raster(system.file("external/test.grd", package="raster"))
plot(r)
 
x<-sampleRandom(r, size=1000, asRaster=TRUE)
plot(x)
 
y<-sampleRandom(r, size=100, sp=TRUE)
plot(r)
plot(y, add=TRUE, pch=20, col="black")
13 Jan

Plot points with a google map background

Good Morning Followers!

Yesterday, I found a nice post here by Joseph Rickert. I produces nice plots with points and a google map background. Try it and have fun!

library(ggmap)
library(plyr)
library(gridExtra)
temp <- tempfile()
download.file("http://www.plantsciences.ucdavis.edu/plant/data.zip", temp)
connection <- unz(temp, "Data/Set3/Set3data.csv")
rice <- read.csv(connection)
names(rice) <- tolower(names(rice))
# Create a custom soil attribute plot
# @param df Data frame containing data for a field
# @param attribute Soil attribute
# @return Custom soil attribute plot
create_plot <- function(df, attribute) {
  map <- get_map(location = c(median(df$longitude), 
                              median(df$latitude)),
                 maptype = "satellite",
                 source = "google",
                 crop = FALSE,
                 zoom = 15)
  plot <- ggmap(map) + geom_point(aes_string(x = "longitude", 
                                             y = "latitude",
                                             color = attribute),
                                  size = 5,
                                  data = df)
  plot <- plot + ggtitle(paste("Farmer", df$farmer, "/ Field", df$field))
  plot <- plot + scale_color_gradient(low = "darkorange", high = "darkorchid4")
  return(plot)
}
ph_plot <- dlply(rice, "field", create_plot, attribute = "ph")
ph_plots <- do.call(arrangeGrob, ph_plot)

6a010534b1db25970b01bb07d2c8c1970d-800wi

13 Jan

Plot a google map satellite image from a place name or address

Good morning Gis addict,
Need a HR satellite image from Google Map for a figure?
Here is an easy way!

 
library(dismo)
 
# enter the address here
x <- geocode('Croix du Sud 2,1348 Louvain‐la‐Neuve')
e <- extent(as.numeric(x[5:8])+ c(-0.001, 0.001, -0.001, 0.001))
g <- gmap(e, type = "satellite")
plot(g)

LLN

12 Jan

RGB colour circle legend

Hello world!
Have you ever needed to make a slick legend for an RGB map or plot such as this:


dum.map.GPPintercompMap

Here is how to do it vectorially in R, and then adding it to your plotRGB map:

add.3C.lgd<- function(x0, y0, scl){
 
  # x0  y0 indicate the position of where to place your legend, scl is used to change its size...
 
  # primary colours
  smp=seq(0,2*pi,pi/180)
  c.x=cos(smp)
  c.y=sin(smp)
  c.x[which(smppi*5/3)]=NA
  c.y[which(smppi*5/3)]=NA
  c.x[rev(which(smp>=pi*5/3))]=cos(smp[which(smp>=pi & smp=pi*5/3))]=sin(smp[which(smp>=pi & smp=pi*1/3 & smp=pi  smp=pi*1/3  smp=pi  smp<=pi*4/3)])+sqrt(3)/2
 
  c.x.g=c.x
  c.y.g=c.y
  c.x.b <- c.x*cos(2/3*pi)-c.y*sin(2/3*pi)+1
  c.y.b <- c.x*sin(2/3*pi)+c.y*cos(2/3*pi)
  c.x.r <- c.x*cos(4/3*pi)-c.y*sin(4/3*pi)+0.5
  c.y.r <- c.x*sin(4/3*pi)+c.y*cos(4/3*pi)+sqrt(3)/2
 
  # secondary colours
  smp1=seq(1/3*pi,2/3*pi,pi/180)
  smp2=seq(pi,4/3*pi,pi/180)
  smp3=rev(seq(2/3*pi,pi,pi/180))
 
  ci.x=c(cos(smp1),cos(smp2)+0.5,cos(smp3)+1)
  ci.y=c(sin(smp1),sin(smp2)+sqrt(3)/2,sin(smp3))
 
  c.x.y=ci.x
  c.y.y=ci.y
  c.x.c <- ci.x*cos(2/3*pi)-ci.y*sin(2/3*pi)+1
  c.y.c <- ci.x*sin(2/3*pi)+ci.y*cos(2/3*pi)
  c.x.m <- ci.x*cos(4/3*pi)-ci.y*sin(4/3*pi)+0.5
  c.y.m <- ci.x*sin(4/3*pi)+ci.y*cos(4/3*pi)+sqrt(3)/2
 
  # draw the individual polygons
  polygon(x = scl*c.x.r +x0,y=scl*c.y.r + y0,col = 'red')
  polygon(x = scl*c.x.b +x0,y=scl*c.y.b + y0,col = 'blue')
  polygon(x = scl*c.x.g +x0,y=scl*c.y.g + y0,col = 'green')
  polygon(x = scl*c.x.y +x0,y=scl*c.y.y + y0,col = 'yellow')
  polygon(x = scl*c.x.m +x0,y=scl*c.y.m + y0,col = 'magenta')
  polygon(x = scl*c.x.c +x0,y=scl*c.y.c + y0,col = 'cyan')
 
  # put some labels
  text(x=x0 + scl*0.5,y=y0 + 2.3*scl,labels = 'Var 1')
  text(x=x0 + scl*2.5,y=y0 - 0.5*scl,labels = 'Var 2',srt=60)
  text(x=x0 - scl*1.5,y=y0 - 0.5*scl,labels = 'Var 3',srt=-60)
}
07 Jan

Dreams come true: a function for RGB plots with GGPLOT2

In love with the ggplot2 look? Sick of using plotRGB to display your RGB raster?
Here comes the rggplot function!

Have a look:

 
rggbplot <- function(inRGBRst,npix=NA,scale = 'lin'){
 
  rgblinstretch <- function(rgbDf){
    maxList <- apply(rgbDf,2,max)
    minList <- apply(rgbDf,2,min)
    temp<-rgbDf
    for(i in c(1:3)){
      temp[,i] <- (temp[,i]-minList[i])/(maxList[i]-minList[i])
    }
    return(temp)
  }
 
  rgbeqstretch<-function(rgbDf){
 
    temp<-rgbDf
    for(i in c(1:3)){
      unique <- na.omit(temp[,i])
      if (length(unique>0)){
        ecdf<-ecdf(unique)
        temp[,i] <- apply(temp[,i,drop=FALSE],2,FUN=function(x) ecdf(x))
      }
    }
    return(temp)
  }
 
  if(is.na(npix)){
    if(ncell(inRGBRst)>5000){
      npix <- 5000
    }
    else{
      npix <- ncell(inRGBRst)
    }
  }
  x <- sampleRegular(inRGBRst, size=npix, asRaster = TRUE)
  dat <- as.data.frame(x, xy=TRUE)
  colnames(dat)[3:5]<-c('r','g','b')
 
  if(scale=='lin'){
    dat[,3:5]<- rgblinstretch(dat[,3:5])
  } else if(scale=='stretch'){
    dat[,3:5]<- rgbeqstretch(dat[,3:5])
  }
 
  p <- ggplot()+ geom_tile(data=dat, aes(x=x, y=y, fill=rgb(r,g,b))) + scale_fill_identity()
 
}
 
b <- (brick(system.file("external/rlogo.grd", package="raster")))*6-600
print(rggbplot(b))

rgb_ggplot2

05 Jan

Plot 3D with plotly

Happy New Gis !

Plotly is a package that easily allows to get interactive, publication-quality plots in your web browser. There are many different plots possible, explore the gallery here to have an idea.

# You can reproduce this figure in R with the following code!
 
# Learn about API authentication here: plot.ly/r/getting-started
 
 
#install packages
install.packages("devtools")
library("devtools")
 
install_github("ropensci/plotly")
 
#upgrade plotly
#devtools::install_github("ropensci/plotly")
 
#Define data
x <- y <- seq(-5, 5, len = 200)
X <- expand.grid(x = x, y = y)
X <- transform(X, z = dnorm(x, -2.5)*dnorm(y) - dnorm(x, 2.5)*dnorm(y))
z <- matrix(X$z, nrow = 200)
 
#Plot
library(plotly)
 
# Find your api_key here: plot.ly/settings/api
py <- plotly(username='username', key='password')
data <- list(
  list(
    z = z, 
    x = x, 
    y = y, 
    type = "surface"
  )
)
layout <- list(
  scene = list(
    xaxis = list(
      title="x",
      gridcolor = "rgb(255, 255, 255)", 
      zerolinecolor = "rgb(255, 255, 255)", 
      showbackground = TRUE, 
      backgroundcolor = "rgb(204, 204, 204)"
    ), 
    yaxis = list(
      title="y",
      gridcolor = "rgb(255, 255, 255)", 
      zerolinecolor = "rgb(255, 255, 255)", 
      showbackground = TRUE, 
      backgroundcolor = "rgb(204, 204, 204)"
    ), 
    zaxis = list(
      title="z",
      gridcolor = "rgb(255, 255, 255)", 
      zerolinecolor = "rgb(255, 255, 255)", 
      showbackground = TRUE, 
      backgroundcolor = "rgb(204, 204, 204)"
    )
  )
)
response <- py$plotly(data, kwargs=list(layout=layout))
url <- response$url
25 Dec

Charting time-series as calendar heat maps

Merry GISmas dearest followers!
A couple of weeks ago, I made a post on heatmaps. Today I came across a nice answer on stackoverflow, it goes like this:

# LOAD THE DATA
stock <- "MSFT"
start.date <- "2006-01-12"
end.date <- Sys.Date()
quote <- paste("http://ichart.finance.yahoo.com/table.csv?s=",
               stock, "&a=", substr(start.date,6,7),
               "&b=", substr(start.date, 9, 10),
               "&c=", substr(start.date, 1,4), 
               "&d=", substr(end.date,6,7),
               "&e=", substr(end.date, 9, 10),
               "&f=", substr(end.date, 1,4),
               "&g=d&ignore=.csv", sep="")    
stock.data <- read.csv(quote, as.is=TRUE)
stock.data <- transform(stock.data,
                        week = as.POSIXlt(Date)$yday %/% 7 + 1,
                        wday = as.POSIXlt(Date)$wday,
                        year = as.POSIXlt(Date)$year + 1900)

The data has the following structure:

> head(stock.data)
        Date  Open  High   Low Close   Volume Adj.Close week wday year
1 2014-12-24 48.64 48.64 48.08 48.14 11437800     48.14   52    3 2014
2 2014-12-23 48.37 48.80 48.13 48.45 23648100     48.45   51    2 2014
3 2014-12-22 47.78 48.12 47.71 47.98 26566000     47.98   51    1 2014
4 2014-12-19 47.63 48.10 47.17 47.66 64551200     47.66   51    5 2014
5 2014-12-18 46.58 47.52 46.34 47.52 40105600     47.52   51    4 2014
6 2014-12-17 45.05 45.95 44.90 45.74 34970900     45.74   51    3 2014

And now this is how we chart it:

library(ggplot2)
ggplot(stock.data, aes(week, wday, fill = Adj.Close)) + 
  geom_tile(colour = "white") + 
  scale_fill_gradientn('MSFT \n Adjusted Close',colours = c("#D61818","#FFAE63","#FFFFBD","#B5E384")) + 
  facet_wrap(~ year, ncol = 1)

Here is how it looks like :-)


Calendar heatmap


16 Dec

Copy to clipboard from R on Ubuntu

Hi young padaGIS,

Install xclip:

sudo apt-get install xclip

Then use this function:

#' @title Copy an object in the clipboard
#' 
#' @description Copy an object in the clipboard}
#' 
#' @param sep The object to copy.
#' @param sep A character to be used as separator for each column of the object
#' @param row.names  Copy row names (default is FALSE)
#' @param col.names Copy column names (default is TRUE)
#' @return copy the object as character in the clipboard
#' @author freecube source:http://stackoverflow.com/questions/10959521/how-to-write-to-clipboard-on-ubuntu-linux-in-r
 
clipboard <- function(x, sep="\t", row.names=FALSE, col.names=TRUE){
     con <- pipe("xclip -selection clipboard -i", open="w")
     write.table(x, con, sep=sep, row.names=row.names, col.names=col.names)
     close(con)
}

Example:

vec <- c(1,2,3,4)
 
clipboard(vec)
clipboard(vec, ",", col.names=FALSE)
clipboard(vec, " ", row.names=TRUE)
05 Dec

Plot heat maps for correlations, confusion matrices etc.

Hi folks!

Wanna present your table in a nice way? Have you ever thought of heat maps?
Here is an example on how to use them:

library(ggplot2)
library(plyr) # might be not needed here anyway it is a must-have package I think in R 
library(reshape2) # to "melt" your dataset
library (scales) # it has a "rescale" function which is needed in heatmaps 
library(RColorBrewer) # for convenience of heatmap colors, it reflects your mood sometimes
nba <- read.csv("http://datasets.flowingdata.com/ppg2008.csv")
head(nba)
nba$Name <- with(nba, reorder(Name, PTS))
nba.m <- melt(nba)
nba.m <- ddply(nba.m, .(variable), transform,rescale = rescale(value))
p <- ggplot(nba.m, aes(variable, Name)) + geom_tile(aes(fill = rescale),colour = "white") +  scale_fill_gradient2(low="#006400", mid="#f2f6c3",high="#cd0000",midpoint=0.5)

The nba data.frame has the folowing structure:

> head(nba)
            Name  G  MIN  PTS  FGM  FGA   FGP FTM FTA   FTP X3PM X3PA  X3PP ORB DRB TRB AST STL BLK  TO  PF
1   Dwyane Wade  79 38.6 30.2 10.8 22.0 0.491 7.5 9.8 0.765  1.1  3.5 0.317 1.1 3.9 5.0 7.5 2.2 1.3 3.4 2.3
2  LeBron James  81 37.7 28.4  9.7 19.9 0.489 7.3 9.4 0.780  1.6  4.7 0.344 1.3 6.3 7.6 7.2 1.7 1.1 3.0 1.7
3   Kobe Bryant  82 36.2 26.8  9.8 20.9 0.467 5.9 6.9 0.856  1.4  4.1 0.351 1.1 4.1 5.2 4.9 1.5 0.5 2.6 2.3
4 Dirk Nowitzki  81 37.7 25.9  9.6 20.0 0.479 6.0 6.7 0.890  0.8  2.1 0.359 1.1 7.3 8.4 2.4 0.8 0.8 1.9 2.2
5 Danny Granger  67 36.2 25.8  8.5 19.1 0.447 6.0 6.9 0.878  2.7  6.7 0.404 0.7 4.4 5.1 2.7 1.0 1.4 2.5 3.1
6  Kevin Durant  74 39.0 25.3  8.9 18.8 0.476 6.1 7.1 0.863  1.3  3.1 0.422 1.0 5.5 6.5 2.8 1.3 0.7 3.0 1.8

And here is the output. Nice isn’t it?


heat_map