2

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?

3
  • Given a contour line A, all you know is that all points from that line to any other contour line B, that does not cross a third contour line must be valued between the values of A and B. You'd normally use this property to interpolate a grid of depths over the space. Would that be something useful to you? What you are trying to do seems to be creating a set of polygons of flat values defined by the contour lines, which doesn't actually represent the underlying depths properly..
    – Spacedman
    Jan 13, 2022 at 16:06
  • Thanks for your thoughts @Spacedman. I don't quite understand what you mean by "proper representation". Each depth polygon would represent a depth interval. As far as I understand those intervals would represent the depths within the model "properly" (of course they are not "real" depths always, there is not enough data). If interpolation would solve the issue, then that is what I should try. I'll look into it tomorrow.
    – Mikko
    Jan 13, 2022 at 17:34
  • Interval is the key word here. Of course it does not work when I try to make polygons out of each line. I need to join two lines and make polygons between them.
    – Mikko
    Jan 14, 2022 at 6:37

1 Answer 1

0

You could use an interpolation package (e.g. fields or gstat) to interpolate z values from the depth contours to a raster and then use st_contour from the stars package to extract polygon "isobands".

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.