maps not created correctly with geom_scatterpie

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

and

lat_lon_province.txt

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

Thank you!!

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.