ggplot_stat_density2d plots for ecological distribution

Question

I'm trying to plot ecological distribution of some species of organisms I'm studying over the Arabian/Persian Gulf. Here is a sample of a code I've tried:

Backround layer

library(ggplot2)
library(ggmap)

nc <- get_map("Persian Gulf", zoom = 6, maptype = 'terrain', language = "English")
ncmap <- ggmap(nc,  extent = "device")

Other layers

  ncmap+
    stat_density2d(data=sample.data3, aes(x=long, y=lat, fill=..level.., alpha=..level..),geom="polygon")+
    geom_point(data=sample.data3, aes(x=long, y=lat))+
    geom_point(aes(x =50.626444, y = 26.044472), color="red", size = 4)+
    scale_fill_gradient(low = "green", high = "red") + scale_alpha(range = c(0.00, 0.25), guide = FALSE)

but, I will like to use the stat_density2d to show the distributions of hundreds of species (which are recorded in columns e.g SP1....SPn) over the water body rather than just displaying latitude and longitude.

Also, is it possible to restrict my heat map to just the water body? I'll appreciate any help and recommendations I can get on this please

It would be beneficial to get some of your data sample.data3. `dput(head(sample.data3, 20))` should give enough to play in `ggplot` with — Vedda, Jan 02 '16 at 18:28
@Amstell, here is a dropbox link to a compressed folder that contains a Rstudi project with all the sample data and polygon i'm working with. [link](https://www.dropbox.com/s/5zssgq4kqlykbf0/Persian%20Gulf.rar?dl=0) Thank You — Hammao, Jan 02 '16 at 21:23
@Hammao This may be an old question but your link to your data doesn't work. Do you still need this answered ? Can you include data using `dput` as Amstell requested ? — steveb, Jul 09 '16 at 06:34
@steveb... here https://www.dropbox.com/s/z2yd6n434r8vj07/Persian%20Gulf.rar?dl=0 — Hammao, Jul 10 '16 at 17:36
@Hammao I think I may have found a solution for your requirement to restrict the heatmap only to the water. What I don't quite understand is this part of your question: "I will [...] show the distributions of hundreds of species (which are recorded in columns e.g SP1....SPn) over the water body rather than just displaying latitude and longitude." In the data (`sample.data`) you provide, all the species are measured at the same coordinates. So, what exactly is that you would like to plot? Some more info on the data set would be helpful. — Felix, Jul 21 '16 at 15:34

Felix · Accepted Answer · 2016-07-25T09:50:54.103

My approach to your question is a pragmatic one: simply put the layer of gulf countries over the heatmap distribution. This crops the heatmap accordingly. Note, however, that the heatmap is still calculated as if it weren't cropped. That means the density calculation is not restricted to the water body only, but it is simply cropped visually.

For the sake of reproducibility, the following code assumes that you've unzipped the .rar file provided by @Hammao and execute the code in the resulting Persian Gulf folder.

# get sample data
sample.data <- read.csv("sample.data3.csv")

Now, we need to get the country shapes for the Gulf countries. I use the rworldmap package for this.

# loading country shapes
library(rworldmap) 

# download map of the world
worldmap <- getMap(resolution = "high") # note that for 'resolution="high"' 
                                        # you also need the "rworldxtra" pkg

# extract Persian Gulf countries...
gulf_simpl <- worldmap[worldmap$SOVEREIGNT == "Oman" | 
                         worldmap$SOVEREIGNT == "Qatar"  |
                         worldmap$SOVEREIGNT == "United Arab Emirates" |
                         worldmap$SOVEREIGNT == "Bahrain" |
                         worldmap$SOVEREIGNT == "Saudi Arabia" |
                         worldmap$SOVEREIGNT == "Kuwait" |
                         worldmap$SOVEREIGNT == "Iraq" |
                         worldmap$SOVEREIGNT == "Iran", ]

# ... and fortify the data for plotting in ggplot2
gulf_simpl_fort <- fortify(gulf_simpl)

# Now read data for the Persian Gulf, which we need to get the distances for
# the extension of the map
PG <- readOGR(dsn = ".", "iho")
PG <- readShapePoly("iho.shp")

PG <- fortify(PG)

Now, it's simply a matter of plotting the layers in the correct order.

# generate plot
ggplot(sample.data) + 

  # first we plot the density...
  stat_density_2d(aes(x = long, y = lat, 
                      fill = ..level..),
                  geom="polygon", 
                  alpha = 0.5) +

  # ... then we plot the points
  geom_point(aes(x = long, y = lat)) +

  # gradient options
  scale_fill_gradient(low = "green", high = "red") + 
  scale_alpha(range = c(0.00, 0.25), guide = FALSE) +

  # and now put the shapes of the gulf states on top
  geom_polygon(data = gulf_simpl_fort, 
               aes(x = long, 
                   y = lat, group = group), 
               color = "black", fill = "white", 
               inherit.aes = F) +

  # now, limit the displayed map only to the gulf 
  coord_equal(xlim = c(min(PG_fort$long), max(PG_fort$long)), 
              ylim = c(min(PG_fort$lat), max(PG_fort$lat))) +
  theme_bw()

this is exactly what i wanted to do... however, instead of just lat and long, i want to plot for each variable... (i'll end up with several plots). lat & long is just to spatial locate the variable. i will modify this answer by adding a loop statement so that the "fill = ..level.." of your answer will accommodate the changing SP1---SPn. I hope this makes the initial question clearer. — Hammao, Jul 23 '16 at 19:47
Ok thanks for the clarification. I've updated the answer. To plot lat / long for each species separately you could also try using facets. Since there are so many species, there will be many facets, however, and reading the plot might become difficult. So looping through the species, as you've suggested, and saving separate plots might be the best way to go, after all. — Felix, Jul 25 '16 at 09:53

ggplot_stat_density2d plots for ecological distribution

1 Answers1