I have a large dataset of depth contours downloaded from an open database. At the time of writing, these contours appear to contain the best available continuous depth model for my region of interest. The second-best continuous alternative is the 15-arc second GEBCO grid, which is not of high enough resolution for my application.
I need non-overlapping valid polygons to use the depth data in my application which uses R and the sf package. I manage to do the raster data polygonization but cannot polygonize the depth contours despite multiple days of trying. I need to do the polygonization for all three files in the downloadable dataset combined but will demonstrate the problem using a subset of the dataset extracted from the main dataset with (for full documentation):
dt <- sf::read_sf("Basisdata_0000_Norge_25833_DybdedataKurverGeneraliserte_Shape/no306300l.shp") %>% dplyr::select(DYBDE, geometry) %>% st_crop(., st_bbox(c(xmin = 17, ymin = 70, xmax = 17.1, ymax = 70.1), crs = 4326) %>% st_as_sfc() %>% st_transform(st_crs(tmp)) ) Packages:
library(sf) #> Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE library(tidyverse) Example data:
dt <- structure( list(DYBDE = c(750, 200, 450, 300, 100, 500, 350, 700, 800, 250, 150, 150, 550, 400, 650, 600), geometry = structure(list( structure(c(575964.768344383, 576066.46, 576282.41, 576393.71, 576491.81, 576570.27, 576714.27, 576747.74810873, 7776554.46439178, 7776744.81, 7777178.12, 7777416.12, 7777644.08, 7777851.88, 7778311.99, 7778399.55686226), .Dim = c(8L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(575964.768344383, 576053, 576540.9, 576904.12, 576978, 577116.19, 577627.56, 577783.23, 577966.37, 578290.68, 578596.55, 578644.52, 578822.28, 578943.06, 579494.32, 579836.51, 580002.01, 580064.87, 580145.9120226, 7772344.2899831, 7772375.77, 7772475.19, 7772640.02, 7772691.13, 7772819.71, 7773381.07, 7773525.77, 7773636.38, 7773783.08, 7773963.91, 7774006.15, 7774198.17, 7774295.27, 7774713.37, 7774926.62, 7775052.29, 7775118.31, 7775226.88181071), .Dim = c(19L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(580145.9120226, 579992.52, 579867.69, 579668.81, 579569.86, 579407.8, 579184.03, 578657.23, 578403.74, 578305.44, 578068.49, 577347.46, 577051.18, 576511.14, 576302.08, 575964.768344383, 7777483.42968411, 7777254.69, 7777045.17, 7776677.81, 7776540.81, 7776366.62, 7776157.47, 7775744.3, 7775513.14, 7775434.18, 7775274.92, 7774836.31, 7774673.32, 7774409.31, 7774286.19, 7774062.91132723 ), .Dim = c(16L, 2L), class = c("XY", "LINESTRING", "sfg" )), structure(c(580145.9120226, 579997.3, 579876.81, 579630.79, 579409.81, 578797.02, 578471.29, 577960.75, 577392.66, 576862.69, 576646.67, 576521.94, 576345.54, 576285.64, 576218.64, 576060.23, 575964.768344383, 7776012.12036657, 7775834.83, 7775708.54, 7775475.33, 7775293.54, 7774813.54, 7774595.12, 7774331.77, 7774014.54, 7773816.51, 7773700.87, 7773598.19, 7773377.26, 7773313.29, 7773252.82, 7773143.16, 7773093.03100332), .Dim = c(17L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(580145.9120226, 580085.69, 579843.63, 579499.66, 579300.43, 579130.25, 578966.74, 578664.36, 578537.64, 578403.31, 578272.59, 578159.32, 578072.44, 578013.93, 577977.54, 577928.39, 577900.41, 577867.97, 577798.6, 577725.91, 577572.65, 577462.73, 577424.06, 577400.93, 577393.13, 577398.06, 577413.43, 577437.77, 577470.71, 577485.158730192, 7767530.39019674, 7767547.76, 7767584.07, 7767597.4, 7767594.63, 7767573.33, 7767521.28, 7767381.15, 7767332.23, 7767294.84, 7767273.08, 7767275.16, 7767303.24, 7767355.52, 7767421.88, 7767546.39, 7767589.77, 7767617.82, 7767652.39, 7767669.18, 7767652.59, 7767602.43, 7767566.12, 7767524.74, 7767472.1, 7767410.27, 7767337.32, 7767252.89, 7767159.14, 7767125.17094982 ), .Dim = c(30L, 2L), class = c("XY", "LINESTRING", "sfg" )), structure(c(575964.768344383, 576013.53, 576516.18, 577074.98, 577247.49, 577622.44, 578453.25, 579006.86, 579143.38, 579235.36, 579464.6, 579654.78, 579954.74, 580145.9120226, 7774401.41705481, 7774431.69, 7774682.87, 7774988.1, 7775094.59, 7775352.93, 7775999, 7776409.32, 7776539.83, 7776665.11, 7777116.55, 7777460.03, 7777956.74, 7778241.72862724), .Dim = c(14L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(575964.768344383, 575982.9, 576102.62, 576357.91, 576584.73, 576778.3, 577248.65, 578111.2, 578421.11, 578729.95, 579533.87, 579713.6, 579866.68, 580098.16, 580145.9120226, 7773413.03712998, 7773423.31, 7773508.76, 7773763.26, 7773923.55, 7774022.94, 7774218.97, 7774675.94, 7774852.87, 7775072.57, 7775708.52, 7775867.82, 7776035.84, 7776347.93, 7776426.23016319), .Dim = c(15L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(577550.647525516, 577532.24, 577464.48, 577413.06, 577337.06, 577188.4, 577024.16, 576771.02, 576671.01, 576543.87, 576238.13, 576128.06, 576076.21, 575964.768344383, 7778399.55686226, 7778341.12, 7778015.55, 7777860.7, 7777701.31, 7777465.73, 7777268.39, 7777008.29, 7776873.24, 7776663.53, 7776132.39, 7775979.68, 7775933.34, 7775877.94694177), .Dim = c(14L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(575964.768344383, 576018.99, 576073.61, 576094.16, 576161.75, 576163.859215472, 7777491.78562878, 7777639.97, 7777816.28, 7777913.66, 7778388.58, 7778399.55686226 ), .Dim = c(6L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(575964.768344383, 576175.56, 576262.23, 576403.49, 576564.62, 576786.1, 576829.04, 576959.16, 576997.57, 577108.84, 577561.43, 577955.7, 578269.11, 578605.59, 578853.45, 579094.05, 579501.19, 579747.51, 579823.65, 579950.65, 580056.08, 580145.9120226, 7772725.43161701, 7772792.35, 7772857.19, 7773018.68, 7773177.11, 7773343.45, 7773385.48, 7773577.63, 7773615.79, 7773677.45, 7773810.96, 7774010.16, 7774218.17, 7774411.82, 7774569.34, 7774740.28, 7775052.41, 7775192.61, 7775263.12, 7775434.27, 7775550.13, 7775631.06921637), .Dim = c(22L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(580145.9120226, 580121.84, 580017.2, 579915.39, 579825.93, 579709.3, 579683.12, 579672.17, 579670.81, 579682.12, 579732.08, 579732.98, 579692.08, 579682.66, 579675.04, 579645.39, 579610.57, 579511.8, 579361.28, 579267.04, 579078.71, 579044.78, 579015.79, 578974.26, 578930.39, 578900.42, 578888.47, 578872.26, 578851.74, 578828.07, 578716.51, 578678.88, 578581.43, 578378.53, 578094.35, 577953.12, 577737.23, 577552.76, 577436.11, 577380.28, 577334.66, 577253.85, 577233.84, 577210.14, 577179.12, 577085.4, 576516.4, 576335.64, 576245.06, 576155.09, 576067.72, 575985.96, 575964.768344383, 7771435.65441595, 7771443.28, 7771501.67, 7771579.38, 7771661.54, 7771795.16, 7771848.51, 7771908.27, 7771981.69, 7772185.73, 7772820.74, 7772954.68, 7773387.24, 7773599.88, 7773630.61, 7773689.73, 7773711.78, 7773642.37, 7773585.14, 7773531.39, 7773391.36, 7773359.2, 7773315.44, 7773201.81, 7773004.72, 7772718.84, 7772631.33, 7772557, 7772494.31, 7772443.29, 7772239.88, 7772188.08, 7772111.31, 7771999.68, 7771773.12, 7771636.63, 7771499.82, 7771325.03, 7771193.12, 7771098.27, 7770957.26, 7770609.89, 7770549.09, 7770503.67, 7770467.92, 7770393.21, 7770002.44, 7769914.97, 7769894.64, 7769891.69, 7769906.18, 7769940.39, 7769955.18395155), .Dim = c(53L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(579528.6, 579558.38, 579584.97, 579609.44, 579632.74, 579676.95, 579735.09, 579760.02, 579765.28, 579765.21, 579747.79, 579709.5, 579656.4, 579599.08, 579549.26, 579531.52, 579520.93, 579519.28, 579528.6, 7771200.01, 7771125.39, 7771081.98, 7771063.29, 7771060.37, 7771086.99, 7771149.12, 7771194.31, 7771219.91, 7771246.87, 7771300.52, 7771351.59, 7771386.66, 7771390.27, 7771357.56, 7771328.13, 7771290.88, 7771246.18, 7771200.01), .Dim = c(19L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(575964.768344383, 576761.81, 577108.44, 577224.93, 577300.92, 577449.37, 577587.92, 577999.71, 578427.6, 578653.19, 578751.01, 578906.08, 579011.34, 579083.32, 579240.03, 579309.04, 579454.28, 579507.37, 579520.88, 579521.270992063, 7774816.69166097, 7775140.44, 7775310.65, 7775384.28, 7775464.3, 7775653.38, 7775793.48, 7776136.81, 7776453.48, 7776641.23, 7776741.58, 7776944.31, 7777120.25, 7777279.87, 7777753.92, 7777913.13, 7778210.73, 7778340.58, 7778393.12, 7778399.55686226 ), .Dim = c(20L, 2L), class = c("XY", "LINESTRING", "sfg" )), structure(c(575964.768344383, 576341.93, 576535.23, 577350.99, 578221.24, 578343.94, 578546.85, 578690, 578892.69, 579241.45, 579483.68, 579540.22, 579646.28, 579859.44, 580105.25, 580145.9120226, 7773744.05616317, 7774039.53, 7774160.28, 7774543.82, 7775028.38, 7775102.26, 7775242.73, 7775367.86, 7775566.32, 7775810.71, 7775997.24, 7776049.71, 7776171.46, 7776469.56, 7776857.78, 7776926.43941226), .Dim = c(16L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(575964.768344383, 576011.21, 576224.95, 576414.97, 576598.24, 576696.76, 577128.52, 577399.64, 577555.3, 577697.78, 577906.21, 577981.94, 578055.51, 578135.42, 578264.110185051, 7775507.69835934, 7775524.23, 7775614.63, 7775731.24, 7775898.16, 7776006.19, 7776529.8, 7776846.96, 7777047.76, 7777259.91, 7777602, 7777771.31, 7777963.5, 7778152.18, 7778399.55686226 ), .Dim = c(15L, 2L), class = c("XY", "LINESTRING", "sfg" )), structure(c(578901.15665015, 578899.67, 578846.57, 578636.96, 578525.43, 578429.89, 578365.53, 578243.57, 578056.46, 577840.56, 577533.83, 577362.66, 577215.11, 577007.64, 576845.75, 576740.41, 576606.21, 576437.23, 576019.21, 575964.768344383, 7778399.55686226, 7778393.24, 7778244.19, 7777793.48, 7777516.6, 7777234.39, 7777119.68, 7776964.6, 7776754.63, 7776548.05, 7776297.54, 7776144.58, 7775986.53, 7775741.25, 7775586.78, 7775511.69, 7775433.48, 7775354.43, 7775193.45, 7775174.21144274), .Dim = c(20L, 2L), class = c("XY", "LINESTRING", "sfg"))), n_empty = 0L, crs = structure(list( input = "WGS 84 / UTM zone 33N", wkt = "PROJCRS[\"WGS 84 / UTM zone 33N\",\n BASEGEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]],\n CONVERSION[\"UTM zone 33N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"Degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",15,\n ANGLEUNIT[\"Degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n ID[\"EPSG\",32633]]"), class = "crs"), class = c("sfc_LINESTRING", "sfc"), precision = 0, bbox = structure(c(xmin = 575964.768344383, ymin = 7767125.17094982, xmax = 580145.9120226, ymax = 7778399.55686226 ), class = "bbox"))), row.names = c(NA, 16L), class = c("sf", "tbl_df", "tbl", "data.frame"), sf_column = "geometry", agr = structure(c(DYBDE = NA_integer_), class = "factor", .Label = c("constant", "aggregate", "identity"))) plot(dt) 
I cannot manage to add the boundaries when casting to polygons:
dt2 <- dt %>% st_cast("LINESTRING") %>% mutate(n = mapview::npts(., by_feature = TRUE)) %>% dplyr::filter(n > 4) # because "polygons require at least 4 points" with a larger dataset dt2 %>% sf::st_cast(.,"POLYGON") %>% plot(max.plot = 1) 
Does not work with st_polygonize() either:
dt %>% st_cast("LINESTRING") %>% st_polygonize() %>% plot() 
I should somehow make sf understand that I want to make the polygons to another polygon encompassing the contour lines:
dt %>% st_boundary() %>% filter(!st_is_empty(.)) %>% st_combine() %>% st_convex_hull() %>% plot() plot(dt, add = TRUE) 
dt %>% st_boundary() %>% filter(!st_is_empty(.)) %>% st_combine() %>% st_convex_hull() %>% st_as_sf() %>% rename("geometry" = "x") %>% bind_rows(., dt2) %>% sf::st_cast(.,"POLYGON") %>% plot(max.plot = 1) 
Created on 2022-01-13 by the reprex package (v2.0.1)
Still not better. Any suggestions how to create depth polygons out of these contours?
