0

I have an issue with the geom_scatterpie function.

I am using R version 4.1.1, scatterpie_0.2.1 and ggplot2_3.4.2 to create some maps with pie charts.

Unfortunately, as you can see in the attached figures, my maps are not being created correctly. map_try_1 map_try_2

This is the script I am running

library(ggplot2)
library(scatterpie)

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

lat_lon_province <- read.csv("lat_lon_province.txt", sep ="\t")
pies <- read.csv("pies.txt", sep ="\t")
colnames(pies)[colnames(pies) == "Province"] ="PROVINCE"

pies_ecoregion_lat_lon <- merge(pies, lat_lon_province, by = "PROVINCE")
pies_ecoregion_lat_lon[is.na(pies_ecoregion_lat_lon)] <- 0
colnms=c("Copepoda", "Isopoda", "Nematoda", "Cestoda", "Trematoda", "Monogenea", "Myxozoa")
pies_ecoregion_lat_lon$Total_parasites<-rowSums(pies_ecoregion_lat_lon[,colnms])

n <- nrow(pies_ecoregion_lat_lon)
pies_ecoregion_lat_lon$radius <- 6 * abs(rnorm(n))

lon_min <- 80
lon_max <- 270
lat_min <- -50
lat_max <- 50

# Using map_data()
pacific_world <- map_data("world2")

#world2MapEnv or world2
#This is an alternative version of the world database based on latitudes [0, 360), 
#which then has the Pacific Ocean in the centre of the map.
mapplot_pacific <- ggplot(pacific_world) + 
  geom_map(data = pacific_world, map = pacific_world, 
           aes(x=long, y=lat, map_id=region), col = "white", fill = "gray50") + theme_bw()
mapplot_pacific + 
  coord_map(xlim=c(lon_min, lon_max), ylim=c(lat_min, lat_max)) +
  geom_scatterpie(aes(x=Longitude, y=Latitude, r=radius), 
                  data = pies_ecoregion_lat_lon, cols = colnames(pies_ecoregion_lat_lon[,c(3:9)]), 
                  color=NA, alpha=.8) +
  geom_scatterpie_legend(pies_ecoregion_lat_lon$radius, x=240, y=35)
ggsave("map_try_1.png", width = 10)

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

pacific_map <- fortify(pacific_world, 
                       xlim = c(lon_min, lon_max), 
                       ylim = c(lat_min, lat_max), 
                       fill = T, plot = F)

pacific <- ggplot(pacific_map)

ggplot(pacific_map, aes(x = long, y = lat, group = group)) + geom_polygon() + 
  coord_map(xlim=c(lon_min, lon_max), ylim=c(lat_min, lat_max)) +
  geom_scatterpie(aes(x=Longitude, y=Latitude, group = PROVINCE, r=radius), 
                  data = pies_ecoregion_lat_lon, cols = colnames(pies_ecoregion_lat_lon[,c(3:9)]), 
                  color=NA, alpha=.8) +
  geom_scatterpie_legend(pies_ecoregion_lat_lon$radius, x=240, y=35) + 
  theme_bw()
ggsave("map_try_2.png", width = 10)

My input files are pies.txt

Province Fish Copepoda Isopoda Nematoda Cestoda Monogenea Trematoda Myxozoa
Eastern Coral Triangle 3 24 6 4 2
Java Transitional 14 1
Lord Howe and Norfolk Islands 2 3 2 1
Northeast Australian Shelf 82 89 28 21 55 109 238 19
Northwest Australian Shelf 51 14 7 3 2 26 2
Sahul Shelf 22 46 26 4 28 12 9
South China Sea 563 24 62 2 2
South Kuroshio 21 24
Sunda Shelf 311 29 11 3 1
Tropical Northwestern Pacific 2 31 14 9
Tropical Southwestern Pacific 14 913 154 48 82 220 393 27
Western Coral Triangle 26 605 122
Central Polynesia 1 42 10 6 2
Easter Island 1 1
Hawaii 27 254 225 32 16 66 668 8
Marquesas 213 4 1
Marshall, Gilbert and Ellis Islands 407 96 10 37 16 105 13
Southeast Polynesia 1041 8 36 1 18 21

and

lat_lon_province.txt

PROVINCE Latitude Longitude
Eastern Coral Triangle -9.989942453 153.1554975
Java Transitional -9.398129883 103.2355709
Lord Howe and Norfolk Islands -30.3708992 163.7689972
Northeast Australian Shelf -14.61130754 147.613617
Northwest Australian Shelf -17.00240261 117.3920957
Sahul Shelf -9.727364469 129.9578054
South China Sea 15.50365482 113.5205013
South Kuroshio 23.64450073 127.5490036
Sunda Shelf 4.935839414 107.017955
Tropical Northwestern Pacific 8.891712166 148.0158917
Tropical Southwestern Pacific -20.82816019 177.8582737
Western Coral Triangle -4.840562086 127.5081903
Central Polynesia -9.675919844 -164.9687772
Easter Island -26.84119987 -107.4769974
Hawaii 22.68499947 -163.0189972
Marquesas -9.251319885 -139.6159973
Marshall, Gilbert and Ellis Islands 4.272244259 171.2996466
Southeast Polynesia -20.34943869 -159.067165

Can you let me know if you have any idea why this is happening?

Thank you!!

Christina P
  • 3
  • 4

1 Answers1

0

This is a CRS issue and not a scatterpie issue. Similar errors are mentioned in other GIS related questions like Appropriate map projection for the Pacific Ocean, How to change map projection from pacific to Atlantic centered?, Reprojecting over the Pacific Ocean with EPSG:3832 gives artefacts.

For the pacific ocean maps the CRS most commonly used is the EPSG:3832, but as mentioned here

the scope of the EPSG:3832 is not for global mapping and its area use can be seen at - http://spatialreference.org/ref/epsg/3832/

So in order to avoid this, first transform a global map to EPSG:3832 and then make valid to neglect the inconsistency.

In R with the sf package this is done with functions st_transform and st_make_valid. And then with st_crop, crop the map to the desired extend. Here is the script with the

#!/usr/bin/env Rscript

library(ggplot2)
library(sf)
library(spData) # for the shapefiles of world
library(scatterpie)

##############################################################################
### load user data and transformations ####
lat_lon_province <- read.csv("lat_lon_province.txt", sep ="\t")
pies <- read.csv("pies.txt", sep ="\t")
colnames(pies)[colnames(pies) == "Province"] ="PROVINCE"

pies_ecoregion_lat_lon <- merge(pies, lat_lon_province, by = "PROVINCE")
pies_ecoregion_lat_lon[is.na(pies_ecoregion_lat_lon)] <- 0
colnms=c("Copepoda", "Isopoda", "Nematoda", "Cestoda", "Trematoda", "Monogenea", "Myxozoa")
pies_ecoregion_lat_lon$Total_parasites<-rowSums(pies_ecoregion_lat_lon[,colnms])

## trasform to sf odject and assign CRS
pies_sf <- st_as_sf(pies_ecoregion_lat_lon, coords=c("Longitude","Latitude"),crs = st_crs("WGS84"))

## transform to EPSG:3832 to match the pacific ocean CRS
pies_sf_3832 <- st_transform(pies_sf, "EPSG:3832")

## make data frame again because the geom_scatterpie doesn't accept sf objects
pies_sf_3832$Longitude  <- st_coordinates(pies_sf_3832)[,1]
pies_sf_3832$Latitude  <- st_coordinates(pies_sf_3832)[,2]
pies_sf_3832_df <- st_drop_geometry(pies_sf_3832)
pies_sf_3832_df$radius <- pies_sf_3832_df$Total_parasites*500

### Load map data from spData and transformations

world <- world # use the spData \


## transform to EPSG:3832 to match the pacific ocean CRS
world_pacific <- st_transform(world, "EPSG:3832") 
## its important to make valid because there are some values 
## in the polygons that go to infinity
world_pacific <- st_make_valid(world_pacific)

## The scope of the EPSG:3832 is not for global mapping 
## therefore keep only the desired area
world_pacific_crop <- st_crop(world_pacific, xmin=-10387065, ymin=-11225295, xmax=19953635, ymax=9906424 )


## plot
##
map_world_pacific_sf <- ggplot() +
    geom_sf(world_pacific_crop, mapping=aes()) +
    geom_sf(pies_sf_3832,  mapping=aes()) +
    geom_label(data = pies_sf_3832_df, 
               mapping=aes(x = Longitude, y = Latitude, label = PROVINCE),
               size = 1.5,
               nudge_x = 1000000,
               nudge_y = 100000)+
    geom_scatterpie(aes(x=Longitude, y=Latitude, r=radius), 
                  data = pies_sf_3832_df,
                  cols = colnames(pies_sf_3832_df[,c(3:9)]),
                  color=NA, alpha=.8) + 
    geom_scatterpie_legend(pies_sf_3832_df$Total_parasites*500,
                           x=-5205792.4,
                           y=-7529715.0) +
    theme_bw()

ggsave("map_pacific_pie_1.png", width = 10, plot=map_world_pacific_sf)

Pacific Ocean map

Note that the aesthetics of radius and legend of geom_scatterpie need some tweaking to get them to the desired form.

Savvas Paragkamian
  • 16
  • 2