ggplot equivalent to locator for maps

[Tyler Rinker]

Hello ggplot community,

I'm using ggplot for creating chloropleths and have the need to plot county names on a state map.  Most names plot beautifully but a few need tweaking.  Tweaking could be made much easier with a locator for ggplot while working with maps.  I put forth this problem at stackoverflow.com and received terrific assistance in solving the problem.  Dwin and Baptise worked together to give me the items to build a function for locator capabilities on a map in ggplot.  It still feels like there's something better as I need to 

[Tyler Rinker]

---

So for the premature send earlier.
==========================

Hello ggplot community,

I'm using ggplot for creating chloropleths and have the need to plot county names on a state map.  Most names plot beautifully but a few need tweaking.  Tweaking could be made much easier with a locator for ggplot while working with maps.  I put forth this problem at stackoverflow.com and received terrific assistance in solving the problem.  Dwin and Baptise worked together to give me the items to build a function for locator capabilities on a map in ggplot.  It still feels like there's something better as it currently only works with one click at a time.  Is there a better way to approach this problem (maybe even something more universal that works outside of maps) or is this it?

Here is a link to the locator code in my drop box as well: http://dl.dropbox.com/u/61803503/gglocator.txt

Thank you in advance.

########################################

require(maps); require(ggplot2); require(grid)  #necessary packages

county_df <- map_data('county')  #mappings of counties by state

ny <- subset(county_df, region=="new york")   #subset just for NYS

ny$county <- ny$subregion

NY <- ggplot(ny, aes(long, lat)) +    #plot a map of NY state with counties

          geom_polygon(aes(group=group), colour='black', fill=NA) +

          coord_map() + geom_point(aes(c(-78, -73), c(41, 40.855), colour=c("blue", "red"))) +

          scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0))

NY 

############FUNCTION_BEGIN#############

ggmap.loc <- function(object){  #ggplot map locator function

   # print(object)

    x <- grid.ls()[[1]][grep("panel-", grid.ls()[[1]])] #locate the panel

    seekViewport(x)

    y <-  as.numeric(grid.locator("npc"))

    locatedX <- min(object$data$long) + y[1]*diff(range(object$data$long))

    locatedy <- min(object$data$lat) + y[2]*diff(range(object$data$lat))

    return(c(locatedX, locatedy))

}

############FUNCTION_END################

ggmap.loc(NY)  #returns coordinates as expected

########################################

[Tyler Rinker]

I had forgot to mention that another down side to this approach is the need to add:

+  scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0))

to the end of your plot.

---

From: tyler_...@hotmail.com
To: ggp...@googlegroups.com
Subject: RE: ggplot equivalent to locator for maps
Date: Mon, 27 Feb 2012 01:14:59 -0500

--
You received this message because you are subscribed to the ggplot2 mailing list.
Please provide a reproducible example: http://gist.github.com/270442
 
To post: email ggp...@googlegroups.com
To unsubscribe: email ggplot2+u...@googlegroups.com
More options: http://groups.google.com/group/ggplot2

[David Kahle]

Hi Tyler -

Nice function!

It didn't work when I copied/pasted directly into R (and I didn't check to see why), but I did some editing and got this perhaps more useful version :

ggmap_loc <- function(object = last_plot(), message = FALSE, xexpand = c(.05,0), yexpand = c(.05, 0)){

  x <- grid.ls(print = message)[[1]]

  x <- x[ grep("panel-", grid.ls(print=message)[[1]]) ] #locate the panel

  seekViewport(x)

  loc <-  as.numeric(grid.locator("npc"))

  xrng <- with(object, range(data[,deparse(mapping$x)]))

  yrng <- with(object, range(data[,deparse(mapping$y)]))    

    

  xrng <- expand_range(range = xrng, mul = xexpand[1], add = xexpand[2])

  yrng <- expand_range(range = yrng, mul = yexpand[1], add = yexpand[2])    

  c(xrng[1] + loc[1]*diff(xrng), yrng[1] + loc[2]*diff(yrng))

}

# df <- data.frame(xvar = 2:10, yvar = 2:10)

# qplot(xvar, yvar, data = df) + annotate(geom = 'point', x = 3, y = 6)

# ggmap_loc()

I'm sure there is  a more efficient way to do this (i.e. pulling from the object scales directly without having to recalculate them), but this works fine.  It can also apply to any ggplot as you see in the minimal example underneath.

Neat!

david.

[Aurelie Cosandey Godin]

Dear all,

I've been struggling with this plot and would appreciate any help.

I am trying to plot a polygon over my geom_points. This is what I've done so far:

    > names(OT1)# my dataset

    [1] "EID"       "latitude"  "longitude" "month"     "year"      "CPUE"      "TSUM" 

    > dim(OT1)

    [1] 2707    7

    > head(OT1)

                    EID latitude longitude month year CPUE TSUM

    1   167-1-1996-1135 67.70000 -61.81667     9 1996    0    0

    2  167-10-1996-1135 67.71667 -59.18333     9 1996    0    0

    3 167-100-1996-1135 67.86667 -59.43333    10 1996    0    0

    4 167-101-1996-1135 67.95000 -59.58333    10 1996    0    0

    5 167-102-1996-1135 68.10000 -59.76667    10 1996    0    0

    6 167-103-1996-1135 67.81667 -59.38333    10 1996    0    0

    OTz<-OT1[with(OT1,OT1$TSUM=="0"),]#selecting only my zeros

    OTc<-OT1[!with(OT1,OT1$TSUM=="0"),]

    #plotting data with ggplot2 (see attached figure)

    v<-ggplot() + geom_point(aes(longitude, latitude, size=TSUM),data= OTc, colour=alpha("red",0.2)) +  facet_wrap(~month, ncol=2)

    v + geom_point(aes(longitude, latitude),size = 1,colour = alpha("black", 0.2), data = OTz) + opts(title="Otter trawl 1996-2011")

I would like to plot the same polygon shape over each of these figures (see polygone shape attached). I've follow the instruction found at URL: http://r.789695.n4.nabble.com/Re-another-question-on-shapefiles-and-geom-point-in-ggplot2-td3424110.html .

 I can plot the polygon, but have a hard time overlaying my geom_points.

    library(rgdal)

    library(ggplot2)

    library(sp)

    library(maptools)

    gpclibPermit()

    div0A <- readOGR(dsn=".", layer="Projections")

    > div0A <- readOGR(dsn=".", layer="Projections")

    OGR data source with driver: ESRI Shapefile 

    Source: ".", layer: "Projections"

    with 1 features and 5 fields

    Feature type: wkbPolygon with 2 dimensions

    > names(div0A);dim(div0A)

    [1] "Id"         "NPAzimutha" "UTM20"      "UTM19"      "AlberEA"   

    [1] 1 5

    > slotNames(div0A) # l

    [1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"

  

    # add the 'id' variable to the shapefile and use it to merge both files

    div0A@data$id = rownames(div0A@data)

    div0A.df <- as.data.frame(div0A)# convert shapefile to dataframe

    div0A_fort <- fortify(div0A,region="id")# fortify to plot with ggplot2 

    head(div0A_fort)

    > head(div0A_fort)

           long      lat order  hole piece group id

    1 -73.50000 78.16666     1 FALSE     1   0.1  0

    2 -73.50000 75.24043     2 FALSE     1   0.1  0

    3 -73.38552 75.04169     3 FALSE     1   0.1  0

    4 -72.95306 74.78239     4 FALSE     1   0.1  0

    5 -70.11000 74.10167     5 FALSE     1   0.1  0

    6 -68.62608 73.72649     6 FALSE     1   0.1  0

    # Merge shapefile and the as.dataframe shapefile

    div0A_merged <- join(div0A_fort,div0A.df, by ="id")

    head(div0A_merged)

    > head(div0A_merged)

           long      lat order  hole piece group id Id NPAzimutha    UTM20  UTM19  AlberEA

    1 -73.50000 78.16666     1 FALSE     1   0.1  0  0     348877 349232.4 349162 348656.4

    2 -73.50000 75.24043     2 FALSE     1   0.1  0  0     348877 349232.4 349162 348656.4

    3 -73.38552 75.04169     3 FALSE     1   0.1  0  0     348877 349232.4 349162 348656.4

    4 -72.95306 74.78239     4 FALSE     1   0.1  0  0     348877 349232.4 349162 348656.4

    5 -70.11000 74.10167     5 FALSE     1   0.1  0  0     348877 349232.4 349162 348656.4

    6 -68.62608 73.72649     6 FALSE     1   0.1  0  0     348877 349232.4 349162 348656.4

    # Plot the shapefile

    ggplot(div0A_merged, aes(long,lat,group=group)) +

      geom_polygon(aes(data=div0A_merged)) +

      geom_path(color="white") + theme_bw()

When I tried this following code as a test, I get an error message: "Error in `[.data.frame`(plot$data, , setdiff(cond, names(df)), drop = FALSE) : 

  undefined columns selected"... 

    p<-ggplot(div0A_merged, aes(long,lat,group=group)) +

      geom_polygon(aes(data=div0A_merged)) +

      geom_path(color="white") + theme_bw()

    

    p + geom_point(aes(longitude, latitude, size=TSUM),data= OTc, colour=alpha("red",0.2)) +  facet_wrap(~month, ncol=2)

THANK YOU!

[Aurelie Cosandey Godin]

Attached my data (OT1.csv) & shapefile (Archives.zip) to help resolve my problem!

Thank you in advance,

[David Kahle]

And actually coming across a use for it in something I am doing (good timing!), this is even better.  Better name, too.

gglocator <- function(n = 1, object = last_plot(), message = FALSE, xexpand = c(.05,0), yexpand = c(.05, 0)){

  

  if(n > 1){

    df <- NULL

    for(k in 1:n){

      df <- rbind(df, gglocator(object = object, message = message, 

        xexpand = xexpand, yexpand = yexpand))

    }

    return(df)

  }

  

  x <- grid.ls(print = message)[[1]]

  x <- x[ grep("panel-", grid.ls(print=message)[[1]]) ] #locate the panel

  seekViewport(x)

  loc <-  as.numeric(grid.locator("npc"))

  xrng <- with(object, range(data[,deparse(mapping$x)]))

  yrng <- with(object, range(data[,deparse(mapping$y)]))    

    

  xrng <- expand_range(range = xrng, mul = xexpand[1], add = xexpand[2])

  yrng <- expand_range(range = yrng, mul = yexpand[1], add = yexpand[2])    

  point <- data.frame(xrng[1] + loc[1]*diff(xrng), yrng[1] + loc[2]*diff(yrng))

  names(point) <- with(object, c(deparse(mapping$x), deparse(mapping$y)))

  point

}

#df <- data.frame(xvar = 2:10, yvar = 2:10)

#qplot(xvar, yvar, data = df) + annotate(geom = 'point', x = 3, y = 6)

#gglocator(4)

david.

On Feb 27, 2012, at 12:38 AM, Tyler Rinker wrote:

[Winston Chang]

It took some work to get all the packages to install and work properly, and the code to run cleanly, but I managed to do it... I'd suggest in the future including code that will all run when copied and pasted into a new R session (this includes loading the data). Also, it helps to know what files go where, and what the working directory should be.

Anyway, I got this to work in ggplot2 version 0.9.0, which should be officially released tomorrow. I think that the problem was solved by the inherit.aes=FALSE argument, which tells the new geom not to inherit the default aesthetic mappings -- those from ggplot(div0A_merged, aes(long, lat, group=group)). The OTc data frame doesn't have those columns, and so when it inherits those mappings and tries to use them, it has problems.

For the new version of ggplot2, I had to change the way color and alpha are specified. There are probably some other small things I forgot...

library(ggplot2)

OT1 <- read.csv('OT1.csv')

OTz<-OT1[OT1$TSUM==0,]#selecting only my zeros

OTc<-OT1[OT1$TSUM!=0,]

# plotting data with ggplot2 (see attached figure)

library(scales)

v <- ggplot(OTc, aes(longitude, latitude, size=TSUM)) +

  geom_point(colour="red", alpha=0.1) + facet_wrap(~month, ncol=2)

v + geom_point(data = OTz,size = 1,colour = "black", alpha=0.2) +

  opts(title="Otter trawl 1996-2011")

library(rgdal)

library(sp)

library(maptools)

gpclibPermit()

div0A <- readOGR(dsn=".", layer="Projections")

names(div0A)

dim(div0A)

library(gpclib)

# add the 'id' variable to the shapefile and use it to merge both files

div0A@data$id = rownames(div0A@data)

div0A.df <- as.data.frame(div0A)# convert shapefile to dataframe

div0A_fort <- fortify(div0A, region="id")# fortify to plot with ggplot2 

head(div0A_fort)

# Merge shapefile and the as.dataframe shapefile

library(plyr)

div0A_merged <- join(div0A_fort, div0A.df, by="id")

head(div0A_merged)

# Graph with the shapefile

p <- ggplot(div0A_merged, aes(long, lat, group=group)) +

  geom_polygon(aes(data=div0A_merged)) +

  geom_path(color="white") + theme_bw()

# With points added

p + geom_point(aes(x=longitude, y=latitude, size=TSUM), data=OTc,

         colour="red", alpha=0.2, inherit.aes=FALSE) + 

  facet_wrap(~month, ncol=2)

On Mon, Feb 27, 2012 at 10:18 AM, Aurelie Cosandey Godin <God...@dal.ca> wrote:

Attached my data (OT1.csv) & shapefile (Archives.zip) to help resolve my problem!

Thank you in advance,

--
You received this message because you are subscribed to the ggplot2 mailing list.
Please provide a reproducible example: http://gist.github.com/270442

To post: email ggp...@googlegroups.com
To unsubscribe: email ggplot2+u...@googlegroups.com
More options: http://groups.google.com/group/ggplot2

[Aurelie Cosandey Godin]

Thank you Chang,

I've just try to download to new version of ggplot2, but it looks like version 0.9.0 is not yet available.

I did run your code, but I still get an error message when I do p+geom_points:

"Error in `[.data.frame`(plot$data, , setdiff(cond, names(df)), drop = FALSE) : 

  undefined columns selected" 

Is this related to the ggplot2 version? Any ideas?

In the meantime, I've also received another suggestion from Felipe Carrillo, which allows me to plot plot the polygon but without unfortunately the facet_wrap(~month, ncol=2) option. 

library(rgdal)
library(ggplot2)
library(sp)
library(maptools)
gpclibPermit() # you won't need this if using the rgeos package

 

shape <- readOGR(dsn=".", layer="Projections")
names(shape);dim(shape)
shape@data$id <- rownames(shape@data)

 #fortify to plot with ggplot2
   shape_fort <- fortify(shape,region="id")

 

  # Import the csv points

  aurelie_points <- read.csv("OT1.csv")
  head(aurelie_points);dim(aurelie_points)
  

OTz<-aurelie_points[with(aurelie_points,aurelie_points$TSUM=="0"),] #selecting only my zeros
    OTc<-aurelie_points[!with(aurelie_points,aurelie_points$TSUM=="0"),]
   Aurelie:

The code could be compacted a little bit but did it this way so you can see what's being done step by step. I will

try to make facet_wrap work when I get a chance. If you want to see the points from OTz on a legend, you could

shorten the limits of the existing legend and make the other legend show below it...just a thought. 

library(gridExtra)
 # I added the points from OTz to this one so you can see both
 # datsets on the same graph...you might also want to add a title to
 # distinguish the different months.
 a <- ggplot(shape_fort, aes(long,lat,group=group)) +
   geom_polygon(fill="lightblue") +
   geom_path(color="green") + geom_point(data=subset(OTc,month==7),
   colour=alpha("red",0.2),aes(longitude,latitude,size=TSUM,group=NULL)) +
    geom_point(data=subset(OTz,month==7),
   colour=alpha("black",0.2),aes(longitude,latitude,size=TSUM,group=NULL)) +
   theme_bw()

 b <- ggplot(shape_fort, aes(long,lat,group=group)) +
   geom_polygon(fill="lightblue") +
   geom_path(color="green") + geom_point(data=subset(OTc,month==8),
   colour=alpha("red",0.2),aes(longitude,latitude,size=TSUM,group=NULL)) +
   geom_point(data=subset(OTz,month==8),
   colour=alpha("black",0.2),aes(longitude,latitude,size=TSUM,group=NULL)) +
   theme_bw()

c <- ggplot(shape_fort, aes(long,lat,group=group)) +
   geom_polygon(fill="lightblue") +
   geom_path(color="green") + geom_point(data=subset(OTc,month==9),
   colour=alpha("red",0.2),aes(longitude,latitude,size=TSUM,group=NULL)) +
   theme_bw()

d <- ggplot(shape_fort, aes(long,lat,group=group)) +
   geom_polygon(fill="lightblue") +
   geom_path(color="green") + geom_point(data=subset(OTc,month==10),
   colour=alpha("red",0.2),aes(longitude,latitude,size=TSUM,group=NULL)) +
   theme_bw()

e <- ggplot(shape_fort, aes(long,lat,group=group)) +
   geom_polygon(fill="lightblue") +
   geom_path(color="green") + geom_point(data=subset(OTc,month==11),
   colour=alpha("red",0.2),aes(longitude,latitude,size=TSUM,group=NULL)) +
   theme_bw()

 f <- ggplot(shape_fort, aes(long,lat,group=group)) +
   geom_polygon(fill="lightblue") +
   geom_path(color="green") + geom_point(data=subset(OTc,month==12),
   colour=alpha("red",0.2),aes(longitude,latitude,size=TSUM,group=NULL)) +
   theme_bw()

 grid.arrange(a,b,c,d,e,f)

 

<OT.png>

[Dennis Murphy]

Hi Aurelie:

ggplot2-0.9.0 is currently the development version of the software,
but as Winston mentioned, it should be sent to CRAN 'real soon now'
and hopefully make it past the security gate. If all goes well, you
should be able to install it tomorrow or the day after from the CRAN
mirrors. If you can't wait that long, then try this, which requires
the devtools and scales packages:

# install.packages('devtools') # uncomment if not installed on your system
library('devtools')
dev_mode() # should reply dev_mode: ON
install_github('ggplot2')
install_github('scales')
library('ggplot2')
library('scales')
<your code here>

dev_mode() basically allows you to install development versions of R
packages not yet released to CRAN and keeps them separate from the
official CRAN versions. When dev_mode is ON, the development version
of packages is used; if OFF, the CRAN versions are used. If dev_mode
is OFF, the install_github() lines will overwrite the old versions,
which you may or may not want.

Dennis

[Aurelie Cosandey Godin]

Thank you Dennis.

[Winston Chang]

Hm, it looks like both the faceting and the inherited group aes were an issue... it was faceting on month, but div0A_merged didn't have a month column, and I guess ggplot2 0.8.9 doesn't like it when you facet on a non-existing variable.

Here's a way to add a month column to the data frame, and replicate each row for each month:

# Get all the months used in OTc

monthdf <- data.frame(month = unique(OTc$month))

# Merge with div0A_merged 

# (replicate each row in div0A_merged for each month)

div0A_merged_month <- merge(div0A_merged, monthdf)

# Graph with the shapefile

ggplot(div0A_merged_month, aes(long, lat, group=group)) +

  geom_polygon() +

  geom_path(color="white") +

  geom_point(data=OTc, aes(x=longitude, y=latitude, size=TSUM),

             colour="red", alpha=0.2, inherit.aes=FALSE) +

  theme_bw() +

  facet_wrap(~ month, ncol=2)

Or of course you could install 0.9.0, as Dennis described.

-Winston

[Aurelie Cosandey Godin]

AAAh! That works like a charm! Thank you very much Winston.

Aurelie Cosandey-Godin
Ph.D. student, Department of Biology
Industrial Graduate Fellow, WWF-Canada

Dalhousie University | Email: god...@dal.ca | Skype: aureliegodinco

----------------------------------------------------------------------------------------------------------------------------------------------------

Want to learn more about sharks in Atlantic Canada? Visit ShARCC! www.atlanticsharks.org