210127_artikkel/gpkg/01_piirkonnale_ruudustike_lisamine.R

#' Piirkonna ruudustike leidmine
#'
#' Etteantud piirkonna geomeetrilise piirjoone ('piir') järele leitakse selle piirkonna EPSG:3301 projektsioonile vastav piirikast ('boundarybox_3301'). Piirikasti järele leitakse selle kastiga kaetud 5x5 km ('epk10t_grid'), 1x1 km ('epk2t_grid') ja 100x100 meetrit ('epk02t_grid') ruudustikud. Samuti leitakse piirkonna piiri sees olevad ja statistilist ning muud infot sisaldavad 5x5 km ('epk10t') ja 1x1 km ('epk2t') sisaldavad ruudu. Andmed salvestatakse GPKG faili kihtidena. Faili nimeks on objekti nimi.
#'
#' @param obj str Objekti nimi.
#' @param pk Piirkonna geomeetriline joon.
#' @param gpkg_home path Salvestatavate GPKG faili asukoht.
#' @param obj_only TRUE/FALSE Kas geomeetriad leitakse ainult objekti või kogu piirikastiga määratud ala jaoks.
#'
#'
piirkonnale_ruudustike_lisamine <- function(obj = NULL, pk = NULL, gpkg_home = "/tmp", obj_only = TRUE) {
  if (is.null(pk) || !sf::st_is_valid(pk)) {
    cat("\nPalun kontrolli geomeetrilise kujundi õigsust.\n")
    return()
  }
  if (is.null(obj)) {
    cat("\nPuudu on objekti nimi.\n")
    return()
  }
  cat(sprintf("\nAlgparameetrid: objekti nimi %s ja GPKG faili kataloog %s\n", obj, gpkg_home))
  postgres <- sprintf(
    "postgres://dbname=\'%s\' host=%s port=%s user=\'%s\' sslmode=%s password=\'%s\' key=\'fid\' srid=4326 type=Polygon checkPrimaryKeyUnicity=\'1\'",
    conf$dbname, conf$host, conf$port, conf$user, conf$sslmode, conf$password
  )
  dsn <- sprintf("%s/%s.gpkg", gpkg_home, obj)
  input_layer_name <- "piir"
  input_layer <- sprintf("%s|layername=%s", dsn, input_layer_name)
  output_layer_name <- "boundarybox_3301"
  tmp_gpkg_file <- tempfile(fileext = ".gpkg")
  # write to gpkg
  sf::write_sf(pk,
    dsn = dsn,
    layer = "piir", driver = "gpkg", append = FALSE
  )

  ## ------------ 2. piirkonna 3301 projektsiooniga piirikast --------------
  ## 2.1  Esmalt leiame milliset 500x500 meetri ruutudega on objekt kaetud.
  # ruut::qgis_algorithm_search_by_word(str = "grid")
  algorithm <- "native:creategrid"
  # cat(qgisprocess::qgis_show_help(algorithm = algorithm))
  result <- qgisprocess::qgis_run_algorithm(
    algorithm = algorithm,
    CRS = "EPSG:3301",
    # EXTENT = '25.454305528,26.259893095,59.454118579,59.714621582 [EPSG:4326]',
    EXTENT = input_layer,
    HOVERLAY = 0,
    HSPACING = 500,
    TYPE = 2,
    VOVERLAY = 0,
    VSPACING = 500,
    OUTPUT = qgisprocess::qgis_tmp_file(ext = ".gpkg")
    # .quiet = TRUE
  )
  result
  pk_grid <- sf::read_sf(qgisprocess::qgis_output(result, "OUTPUT"))
  # sf::st_geometry(pk_grid) %>% plot()
  ## 2.2 Leitud ruutudele leiame piirikasti.
  algorithm <- "qgis:minimumboundinggeometry"
  # cat(qgisprocess::qgis_show_help(algorithm = algorithm))
  result <- qgisprocess::qgis_run_algorithm(
    algorithm = algorithm,
    FIELD = "",
    INPUT = pk_grid,
    TYPE = 3,
    OUTPUT = tmp_gpkg_file
    # .quiet = TRUE
  )
  result
  # pk_grid_bb <- sf::read_sf(qgisprocess::qgis_output(result, "OUTPUT"))
  # sf::st_geometry(pk_grid_bb) %>% plot()
  # sf::st_geometry(pk) %>% plot(add = T)

  system(sprintf("ogr2ogr -f GPKG -overwrite  %s %s -nln %s -t_srs \"EPSG:3301\"", dsn, tmp_gpkg_file, output_layer_name))

  ## ------------- 3. piirkonna epk10t (5x5 km) ruudud --------------------
  # 3.1 kogu ruutvõrgustik
  ruudud <- c("epk10t_grid", "epk2t_grid", "epk02t_grid")
  j <- 4
  for (j in 1:length(ruudud)) {
    ruut <- ruudud[j]
    output_layer_name <- ruut
    conf <- ruut::get_config()
    pg <- ruut::construct_ogr2ogr_PG_connect_str(conf = conf)
    input <- sprintf(
      '%s table=\"%s\".\"%s\" (geometry)',
      postgres, "maaamet", ruut
    )
    # ruut::qgis_algorithm_search_by_word(str = "extract")
    algorithm <- "native:extractbylocation"
    # cat(qgisprocess::qgis_show_help(algorithm = algorithm))
    result <- qgisprocess::qgis_run_algorithm(
      algorithm = algorithm,
      INPUT = input,
      INTERSECT = sprintf("%s|layername=%s", dsn, "boundarybox_3301"),
      OUTPUT = tmp_gpkg_file,
      PREDICATE = c(0, 1)
      # .quiet = TRUE
    )
    result
    # assign(ruut, sf::read_sf(qgisprocess::qgis_output(result, "OUTPUT")))
    system(sprintf("ogr2ogr -f GPKG -overwrite  %s %s -nln %s -t_srs \"EPSG:3301\"", dsn, tmp_gpkg_file, output_layer_name))
  }
  # sf::st_geometry(epk10t_grid) %>% plot(add = T)

  # 3.2 ainult piirkonnaga seotud ning informatsiooni sisaldav ruutvõrgustik
  ruudud <- c("epk10t", "epk2t")
  j <- 4
  for (j in 1:length(ruudud)) {
    ruut <- ruudud[j]
    output_layer_name <- ruut
    conf <- ruut::get_config()
    # pg <- ruut::construct_ogr2ogr_PG_connect_str(conf = conf)
    input <- sprintf(
      '%s table=\"%s\".\"%s\" (geometry)',
      postgres, "maaamet", ruut
    )
    # ruut::qgis_algorithm_search_by_word(str = "extract")
    algorithm <- "native:extractbylocation"
    # cat(qgisprocess::qgis_show_help(algorithm = algorithm))
    result <- qgisprocess::qgis_run_algorithm(
      algorithm = algorithm,
      INPUT = input,
      INTERSECT = input_layer,
      OUTPUT = tmp_gpkg_file,
      PREDICATE = c(0, 1)
      # .quiet = TRUE
    )
    result
    # assign(ruut, sf::read_sf(qgisprocess::qgis_output(result, "OUTPUT")))
    system(sprintf("ogr2ogr -f GPKG -overwrite  %s %s -nln %s -s_srs \"EPSG:4326\" -t_srs \"EPSG:3301\"", dsn, tmp_gpkg_file, output_layer_name))
  }
  # sf::st_geometry(epk10t) %>% plot(add = T)
  ## ---------------------- vaata layer'id ----------------------
  # Vaata layer'eid
  sf::st_layers(dsn = dsn)
}
# piirkonnale_ruudustike_lisamine(obj = NULL, pk = NULL, gpkg_home = "/tmp")