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

Answer 1

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()