0

I'm trying to extract monthly composite sea ice concentration from the regions of the Southern Ocean in this figure:

enter image description here

I've found the data that I want here. The issue is that the coordinates for those data are in polar stereographic projection, while the rest of my R project is in lat-lon.

Is there a simple way to convert all polar stereographic coordinates to lat-lot? Or, is there an R package that can do that? I have the code to access and extract a subset of the data that I want within a defined polygon, but the regions in the above figure are lat-lon defined polygons, so my workflow will be:

  1. Define lat-lon boundaries of the extraction polygon
  2. Convert to polar coordinates
  3. Extract the data
  4. Convert x-y coordinates of the data to lat-lon
  5. Plot on the figure

I've found a few helpful articles but nothing that can show me how to do that specifically, although this article was close. This article was also very close, but I need the actual underlying data, not just making the map. Both of those are in python though and I'm working in R.

Here is the code for the Antarctica figure:

# Defines the x axes required
x_lines <- seq(-120, 180, by = 60)
wm <- map_data("world")

# CCAMLR Areas 88.1, 88.2, 88.3
Area_88.1 <- data.frame(
  long = c(150, 190, 190, 150),
  lat = c(-60, -60, -85, -85))

Area_88.2 <- data.frame(
  long = c(-170, -105, -105, -170),
  lat = c(-60, -60, -85, -85))

Area_88.3 <- data.frame(
  long = c(-105, -70, -70, -105),
  lat = c(-60, -60, -76, -76))

# Whole Map
ggplot() +
  
  # Longitude lines, dashed. Move below polygon code to make these appear on top.
  geom_segment(aes(y = -57, yend = -90, x = x_lines, xend = x_lines), linetype = "dashed") +
  
  # Adding CCAMLR Area 88
  geom_polygon(data = Area_88.1, aes(x = long, y = lat), fill = "blue", colour = "black", alpha = 0.5) +
  geom_polygon(data = Area_88.2, aes(x = long, y = lat), fill = "red", colour = "black", alpha = 0.5) +
  geom_polygon(data = Area_88.3, aes(x = long, y = lat), fill = "green", colour = "black", alpha = 0.5) +
  
  # Adding whole world land polygon
  geom_polygon(data = wm, aes(x = long, y = lat, group = group), fill = "grey", colour = "black", alpha = 1) +
  
  # Lat + Long labels
  geom_text(aes(x = 150, y = seq(-65, -85, by = -10), hjust = -0.2, label = paste0(seq(65, 85, by = 10), "°S"))) +
  geom_text(aes(x = x_lines[-c(2, 5)], y = -56, vjust = -0.5, label = c("120°W", "0°", "60°E", "180°W"))) +
  geom_label(aes(x = x_lines[c(2, 5)], y = -60, vjust = 0, hjust = 0, label = c("60°W", "120°E")), fill = NA, label.size = 0) +
  
  # Dark lines
  geom_hline(aes(yintercept = -57), linewidth = 1) +
  geom_segment(aes(y = -57, yend = -90, x = 120, xend = 120), colour = "black", size = 1) +
  geom_segment(aes(y = -57, yend = -90, x = -60, xend = -60), colour = "black", size = 1) +
  
  # Zooming in on Antarctica, fixing projection. X is lat, Y is long.
  scale_y_continuous(limits = c(-90, -56), breaks = seq(-45, -90, by = -10), labels = NULL, expand = c(0, 0)) +
  scale_x_continuous(breaks = NULL, expand = c(0, 0)) +
  coord_map("ortho", orientation = c(-90, 210, 0), xlim = c(-60, -240)) +
  
  # Style
  theme(
    panel.background = element_blank(),
    panel.grid.major = element_line(
      linewidth = 0.25,
      linetype = "dashed", colour = "black"),
    axis.ticks = element_blank()) +
  labs(x = NULL, y = NULL)
blitz1259
  • 47
  • 4

1 Answers1

1

If you have data that uses a polar stereographic projection and you need to convert it to latitude/longitude, you will need to use coordinate system transformation functions. In R, one of the packages that provides this capability is the rgdal package.

Here's a simple example of how to convert the coordinates:

library(rgdal)

# Define your polar stereographic projection string (EPSG:3031)
polar_proj <- "+init=epsg:3031"

# Define your latitude/longitude projection string (EPSG:4326)
latlon_proj <- "+init=epsg:4326"

# Convert from polar to lat/lon
# Assume df_polar is a data.frame with columns x, y representing polar coordinates
df_polar <- data.frame(x = c(-3944331, 4075000), y = c(-4237000, -3250000)) # example polar coordinates
coordinates(df_polar) <- ~x+y
proj4string(df_polar) <- CRS(polar_proj)

# Transform to lat/lon
df_latlon <- spTransform(df_polar, CRS(latlon_proj))

# Now df_latlon contains the latitude/longitude coordinates
print(df_latlon)

In your case, you need to first define the polygon boundaries in lat-lon, convert those to polar coordinates, extract the data, and then convert the data coordinates to lat-lon for plotting. You can use the above example as a guide, adjusting as needed for your specific project.

megaultron
  • 399
  • 2
  • 15