Search code examples
rgisr-sfr-sptmap

R: Convert sf polygon defined using lat/long to UTM


I have the following polygon, defined using degrees latitude/longitude:

## Define latitude/longitude
lats <- c(64.25086, 64.24937, 63.24105, 63.22868)
lons <- c(-140.9985, -136.9171, -137.0050, -141.0260)
df <- data.frame(lon = lons, lat = lats)


polygon <- df %>% 
  ## EPSG 3578; Yukon Albers projection
  st_as_sf(coords = c('lon', 'lat'), crs = 3578) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast('POLYGON')

When I plot it on a map using Tmap, it appears in the Pacific Ocean off the coast of British Columbia, rather than in the middle of the Yukon:

library(sf)
library(sp)
library(tmap)
library(dplyr)
library(magrittr)
library(leaflet)

m <- tm_shape(data$study_boundary) + tm_borders(col = 'black',
                                                lwd = 5,
                                                zindex = 1000)

m

I am guessing that the problem is in using lat/long rather than UTMs because I have other polygons defined using UTMs that do appear where they (and the polygon defined above) are supposed to be. I found several other posts going the other way (UTM to lat/long) using spTransform, but I haven't been able to go lat/long to UTM with spTransform. I tried the code below:

poly_utm <- st_transform(polygon, crs = "+proj=utm+7")

But that didn't work either.

Thanks!


Solution

  • This (which I've improved by removing the pipe):

    st_as_sf(df, coords = c('lon', 'lat'), crs = 3578)
    

    creates a spatial points data frame using the numbers in the data frame for the coordinates, and the crs code of 3578 as the label for what those numbers represent. It does not change the numbers.

    It looks like those numbers are actually lat-long coordinates, which means they are probable crs code 4326, the lat-long system used for GPS, also known as WGS 84. But it might not be. But probably is. Do check. So anyway, you should do:

    df_unprojected = st_as_sf(df, coords = c('lon', 'lat'), crs = 4326)
    df_projected = st_transform(df_unprojected, 3578)
    

    The st_transform function does the actual change of the coordinate numbers and assigns the new CRS code to the spatial data metadata. That should give you a set of points you can then plot and check they are in the right place before you throw it into summarise and st_cast.