0

I have a question regarding converting spatial data in R and bringing it from R into QGIS.

I have a GeoTiff of Antarctic sea ice concentration, downloaded from the link below:

https://seaice.uni-bremen.de/databrowser/#day=13&month=10&year=2022&img=%7B%22image%22%3A%22image-1%22%2C%22product%22%3A%22AMSR%22%2C%22type%22%3A%22visual%22%2C%22region%22%3A%22Antarctic3125%22%7D

I want to extract the contour of the sea ice edge (defined as 15%), and then have this contour in a file type that I can open in QGIS and reproject for use in making other maps. My current understanding is that to do this, I would need to convert the contour to a spatial points df, and then convert that to a spatial polygons df which I would then be able to open as a shapefile in QGIS. However, I think I'm going wrong here as I cannot make the conversion with the below code - any suggestions?

**This is my current workflow:**

library(raster)
library(tidyverse)
library(sp)
library(sf)

#Load in sea ice geotiff
sic <- raster('Environmental_Data/SIC/AMSR2/asi-AMSR2-s3125-20220107-v5.4.tif')/1 
plot(sic)

#Make all values over land NA
sic[sic>100] = NA

#Crop to make area smaller (I have a specific area of interest)
sic = crop(sic, extent(sic)*c(0.5,0.5,0,1)) 
plot(sic)

#Pull out the sea ice edge (15% contour) (this makes it a spatial lines df)
ie = rasterToContour(sic, levels=15)

#Convert to spatial points
ie.pt = as(ie, "SpatialPointsDataFrame") plot(ie.pt, add=T, pch=16, cex=0.4)

#Convert to spatial polygons
ie.pt_poly <-as(ie.pt, "SpatialPolygons")

#Then I get this error: 
Error in as(ie.pt, "SpatialPolygons"):
no method or default for coercing “SpatialPointsDataFrame” to “SpatialPolygons”
Phil
  • 7,287
  • 3
  • 36
  • 66
Ellie_Petrels
  • 17
  • 6

2 Answers2

1

reworking your process to terra and sf...

library(terra)
library(sf)
sic <- rast('~/Downloads/asi-AMSR2-s3125-20221113-v5.4.tif')
sic[sic>100] = NA
sic2 = crop(sic, ext(sic)*c(0.5,0.5,0,1))
sic2_contour <- terra::contour(sic2, maxcells=100000, filled = TRUE) # plot side effect
sic2_cont <- as.contour(sic2)
sic2_cont_disagg <- disagg(sic2_cont)
y <- sf::st_as_sf(sic2_cont_disagg)
y
Simple feature collection with 6519 features and 1 field
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: -1968608 ymin: 206384.6 xmax: 1968608 ymax: 3660462
Projected CRS: WGS 84 / NSIDC Sea Ice Polar Stereographic South
First 10 features:
   level                       geometry
1     10 LINESTRING (-1968608 340765...
2     10 LINESTRING (-1955825 232458...
3     10 LINESTRING (-1968608 259539...
4     10 LINESTRING (-1968608 262189...
5     10 LINESTRING (-1960827 264530...
6     10 LINESTRING (-1968608 265308...
7     10 LINESTRING (-1968608 278293...
8     10 LINESTRING (-1943042 251270...
9     10 LINESTRING (-1943042 275001...
10    10 LINESTRING (-1930259 331948...
unique(y$level)
 [1]  10  20  30  40  50  60  70  80  90 100
y10 <- y[which(y$level == 10),]
plot(sic2)
plot(y10, col = 'pink', lwd =3, add = TRUE)

I can't think why one would go to points, except perhaps to then buffer and fill one's contour. But terra::writeRaster(sic2..., and terra::writeVector(y,..., or y10 pull into QGIS and see.

enter image description here

There are pink(s) interior to sic2 as these presumably are holes in sea ice that have the same value as northernmost contour that could perhaps be further removed by testing for within.

Chris
  • 1,647
  • 1
  • 18
  • 25
  • Thank you for your answer - However, I'm trying to get the margin as a single line. The terra code in your answer plots all of the 15% contours, but I'm trying to retain only the outermost (most northerly) 15% delineation. – Ellie_Petrels Nov 14 '22 at 18:11
  • Lengthy discussion on this and short [terra solutions](https://stackoverflow.com/questions/67579831/rasterrastertocontour-contour-lines-are-not-continuous), that I'll use above, but appears to address both your `raster` and `terra`. – Chris Nov 14 '22 at 18:54
0

I think this is what you are looking for.

library(terra)
r <- rast("asi-AMSR2-s3125-20221113-v5.4.tif")

# crop to the area of interest
e <- ext(-1975000, 1975000, 2e+05, 4350000)
re <- crop(r, e)

# get contour and save to file
v <- as.contour(re, levels=15)
writeVector(v, "contour_lines.shp")

Contours are normally lines (neither points nor polygons). But if you wanted a polygon you could do

x <- ifel(x <15 | x>100, NA, 1)
p <- as.polygons(x)
writeVector(p, "contour_polygons.shp")

Or, more generally, use terra::classify to create regions before using as.polygons.

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • This worked, with the inclusion of some lines to plot the .shp in ggplot: ``` R <- terra::rast('asi-AMSR2-s3125-20220107- v5.4.tif')/1 #crop to the area of interest e <- ext(-1975000, 1975000, 2e+05, 4350000) R_Crop <- crop(tp_max, e) #get 15% contour and save to file R_v <- as.contour(R_crop, levels=15) writeVector(R_v, "15_contour_lines.shp") #prep for geom_sf R_sf <- R_v %>% st_as_sf() R_sf <- R_sf %>% st_transform(target_crs) #plot to check ggplot()+ geom_sf(aes(), data = tp_max_sf) ``` – Ellie_Petrels Nov 19 '22 at 17:48