0

I would like to plot shapefiles from the EUFORGEN website (http://www.euforgen.org/) in R (European beech and Scots pine). At first I had problems, because they are files with holes. I was able to solve this problem for beech using sebdalgarno’s hint.

For pine it is a bit more difficult, because with the same method, suddenly only single fragments (isolated population) are displayed but not the core population (see picture).

I could imagine that I could "adjust" something when reading in the file, so that the whole layer is displayed.

Any ideas?

Edit (referring to the comment from @Spacedman)

I'm sorry, but I don't know how to read a shape-file via r-script from the internet. If someone knows how to do that, I can change that with pleasure. Therefore here again the link to the cloud

Cloud_Beech(Fagus sylvatica) and Pine(Pinus sylvestris)

or to the zip-file from the Webpage EUFORGEN.org

Here is the map as it is almost correct. The darker small areas (e.g., in Norway and Sweden) are holes in the polygon, but are not displayed correctly with geom_polygon().

library("rnaturalearth")
library("rgdal")
library("ggplot2")


gg_world <- ne_countries(scale = "medium", returnclass = "sp")  # world map as sp class

Pinus <- readOGR("./chorological_maps_dataset/Pinus sylvestris/shapefiles/Pinus_sylvestris_plg_clip.shp")


ggplot()+
  geom_polygon(data = gg_world, mapping=aes(x = long, y = lat, group = group), fill = "white", color = "grey50")+
  # Species shape file
  geom_polygon(data = Pinus, aes(x = long, y = lat, group = group), colour = "black", fill = "darkblue", alpha = 0.2)+
  coord_cartesian(xlim=c(-15, 40), ylim=c(35, 75), default = F)+
  # Layout
  theme_bw()+
  # axes
  scale_x_continuous(expand = c(0, 0), name = "Longitude (°E)")+
  scale_y_continuous(expand = c(0, 0), name = "Latitude (°N)")

Pine distribution with geom_polygon()

So I continued with sebdalgarno's suggestion and used geom_sf() and as also above, with the hint from Spacedman to use a world map in sf format. I used ne_countries from the rnaturalearth packages

gg_world <- ne_countries(scale = "medium", returnclass = "sf")
Pinus_st <- st_read("./chorological_maps_dataset/Pinus sylvestris/shapefiles/Pinus_sylvestris_plg_clip.shp")


ggplot() + 
 
  
  geom_sf(data = Pinus_st, colour = NA, fill = "darkblue", alpha = 0.5)+
  #geom_sf(data = Fagus_st, colour = NA, fill = "red", alpha = 0.5)+
  geom_sf(data = gg_world, fill = NA, color = "grey50")+
  
  coord_sf(xlim=c(-15, 40), ylim=c(35, 75), expand = F)+
  
  # Layout
  theme_bw()+
  
  # axes
  scale_x_continuous(expand = c(0, 0), name = "Longitude (°E)")+
  scale_y_continuous(expand = c(0, 0), name = "Latitude (°N)")

and get this

only single polygones

A new step-by-step build of the code then showed that the map section boundary was causing the error. As was also the case here Link. If one takes out the alpha command, it works.

gg_world <- ne_countries(scale = "medium", returnclass = "sf")
Pinus_st <- st_read("./chorological_maps_dataset/Pinus sylvestris/shapefiles/Pinus_sylvestris_plg_clip.shp")


ggplot() + 
 
  
  geom_sf(data = Pinus_st, colour = NA, fill = "darkblue")+
  geom_sf(data = gg_world, fill = NA, color = "grey50")+
  coord_sf(xlim=c(-15, 40), ylim=c(35, 75), expand = F)+
  # Layout
  theme_bw()+
  # axes
  scale_x_continuous(expand = c(0, 0), name = "Longitude (°E)")+
  scale_y_continuous(expand = c(0, 0), name = "Latitude (°N)")

Zoom in distribution pine

3
  • Can you show us what you expect the map to look like? Have you tried plotting your data with a simple base "plot" call to see if its a problem with ggplot's maps or with the data? Can you tell us exactly which of the shapefiles in the linked downloads you are working with in order for us to replicate your workflow? Can you explain where getMap comes from because its not in your library list? Can you give Pinus_st <- st_read() an argument so it actually runs? Please edit your question down to a minimal example that runs and illustrates the problem.
    – Spacedman
    Commented Jan 24, 2022 at 14:46
  • Now gg_world appears in your code but its not created anywhere. Also, two of your shapefile links need an OwnCloud password, and the other is to just the .shp file, which won't work properly without the .shx, .prj and .dbf. All you need to say is "Download and unzip chorological_maps_dataset_20211007.zip (Mb) from [link]" and then use paths relative to that, eg Pinus = st_read("./chorological_maps_dataset/Pinus sylvestris/shapefiles/Pinus_sylvestris_plg_clip.shp")
    – Spacedman
    Commented Jan 25, 2022 at 9:18
  • It looks like either the shapefile or ggplot's geom_polygon is broken. I'd switch to sf spatial objects and then use tmap for mapping, or if you want to stick with ggplot, then use geom_sf - but you'll probably have to find a sf country outline basemap because the one from rworldmap is sp class... that could be converted.. except it has some unclosed polygon errors...
    – Spacedman
    Commented Jan 25, 2022 at 9:20

1 Answer 1

0

Something is buggy somewhere. Here's how I'd map that shapefile with tmap:

library(sf)
library(tmap)
library(tmaptools)

data(World)

Pinus = st_read("./chorological_maps_dataset/Pinus sylvestris/shapefiles/Pinus_sylvestris_plg_clip.shp")

xlim=c(-15, 40)
ylim=c(35, 75)


map = tm_shape(World, bbox=bb(xlim=xlim, ylim=ylim)) +
    tm_graticules() + 
    tm_polygons("white") +
    tm_shape(Pinus) + tm_polygons("darkblue", alpha=0.2) 

print(map) # or just "map" if interactive

Adjust colours, labels etc as per requirements, but this is close enough:

enter image description here

I see the holes in Scandinavia which I guess were the real problems.

Or with ggplot, using the same Pinus read as before, as an sf object:

ggplot() + geom_sf(data=World) + geom_sf(data = Pinus, fill="darkblue", alpha=0.2) +
    coord_sf(xlim=xlim, ylim=ylim)

I prefer tmap because its easy to change the mode and make an interactive browser map, but I give the ggplot map as well anyway.

enter image description here

And for completion here's the base graphics version:

plot(st_geometry(World), xlim=xlim, ylim=ylim)
plot(Pinus, add=TRUE, col="#00008040")

enter image description here

1
  • I like your solution with tmap. Thanks for your help!
    – Baloo
    Commented Jan 25, 2022 at 11:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.