ruut/R/pk_lisame_rahvaarvud.R

#' Piirkonnale rahvastikuandmete lisamine
#'
#' Statistikaameti 1x1 km olevate rahvastikuandmete lisamine piirkonna kaardiruutudele. Andmed salvestatakse postgis andmebaasi.
#'
#' @param obj str Objekti nimi. Edaspidi on oluline ainult see nimi. Piirkonna geomeetrilist joont ei ole vaja lisada.
#' @param conf A list() of configuration variables. Default values \code{\link[ruut]{get_config}}.
#' @return Rahvastikuruutude loomine.
#' @seealso [sf::st_read()], [sf::write_sf()],[sf::st_transform()],[ruut::pk_sellest_alustame_db_loomist()],[ruut::pk_lisame_ruudustikud()] ,[ruut::pk_lisame_rahvaarvud()],[ruut::pk_lisame_jooned()],[ruut::pk_lisame_punktid()],[ruut::pk_teisendame_polygoone()],[ruut::pk_teisendame_jooni()],[ruut::pk_teisendame_punkte()]
#' @keywords postgis, boundary box, EPSG:3301
#' @export
#' @examples
#' \dontrun{
#'
#' obj <- "marja"
#' pk_lisame_rahvaarvud(obj = obj, conf = NULL)
#'
#' # Layers list.
#' ruut::db_schema_tablenames(conf = conf)
#' }
pk_lisame_rahvaarvud <- function(obj = NULL, conf = NULL) {
  ## ------------- muutujad ja teisendused ---------------
  vars <- ajutised_muutujad(pk = NULL, obj, conf)
  obj <- vars$obj
  piir <- vars$pk
  conf <- vars$conf
  ## Konfiguratsiooni muutujale väärtuste omistamine, kui seda pole antud.
  if (!any("a00_piir" %in% ruut::db_schema_tablenames(conf = conf))) {
    cat(
      "\nAndmebaas loomata. Palun funktsiooniga ruut::pk_sellest_alustame_db_loomist() andmebaasi loomist.\n"
    )
    return(NULL)
  }
  ## Tabelid mis on juba andmebaasis ja mida ümber ei kirjutata.
  olemasolevad_tabelid <- ruut::db_schema_tablenames(conf = conf)

  ## -------------------- Loop -----------------------
  intersect_layers <- c("bb2_epk2t_grid")
  andmed <-
    data.frame("schema" = character(0), "table" = character(0))
  andmed <-
    rbind(andmed,
          data.frame("schema" = "statistikaamet", "table" = "rel_1x1km"))

  for (intersect in intersect_layers) {
    conf$table <- intersect
    conf$schema <- obj
    input <-
      ruut::construct_to_gpkg_output_postgres_str(
        conf = conf,
        geometry_type = "MultiPolygon",
        srid = 3301,
        checkPrimaryKeyUnicity = TRUE,
        key = "id"
      )

    for (i in 1:nrow(andmed)) {
      ## Kui olemasolev tabel eksisteerib andmebaasis, siis jätame arvutused
      ## selle tabeliga vahele.
      uus_tabel <-
        sprintf("tif_stat_%s_epk02t_tif", andmed$table[i])
      if (uus_tabel %in% olemasolevad_tabelid) {
        cat(sprintf(
          "\nTabel on %s.%s juba andmebaasis olemas.\n",
          obj,
          uus_tabel
        ))
        next
      }
      geom <- "geom"
      # if (andmed$schema[i] %in% c("teeregister_wfs"))  geom <- "geometry" else geom <- "geom"
      ## 1.1 Ruutude sidumine asukoha järele.
      # ruut::qgis_algorithm_search_by_word("fix")
      algorithm <- "native:joinattributesbylocation"
      conf$table <- andmed$table[i]
      conf$schema <- andmed$schema[i]
      intersect_layer <-
        ruut::construct_to_gpkg_output_postgres_str(
          conf = conf,
          geometry_type = "Polygon",
          srid = 3301,
          checkPrimaryKeyUnicity = TRUE,
          key = "id",
          geometry_field = geom
        )
      output <- vars$tmp_gpkg_file_output_1 # ajutine fail
      str <- sprintf(
        "{    'DISCARD_NONMATCHING' : true,
    'INPUT' : '%s',
    'JOIN' : '%s',
    'JOIN_FIELDS' : ['rahvaarv'],
    'METHOD' : 1,
    'OUTPUT' : '%s',
    'PREDICATE' : [1],
    'PREFIX' : ''}",
        input,
        intersect_layer,
        output
      )
      cmd <-
        ruut::construct_qgis_output_result_to_better_format(str = str, algorithm = algorithm)
      system(cmd)
      ## 1.2 Asendame NULL väärtustega 0.
      algorithm <- "qgis:advancedpythonfieldcalculator"
      input <- vars$tmp_gpkg_file_input_1 # ajutine fail
      output <- vars$tmp_gpkg_file_output_2
      cmd <-
        sprintf(
          "qgis_process run %s --FIELD_NAME='value' --FIELD_TYPE=1 --FIELD_LENGTH=8 --FIELD_PRECISION=2 --GLOBAL=\"def getValue(x):
    if not x:
        value = 0
    else:
        value = round(x,1)
    return value\" --FORMULA='value = getValue( <rahvaarv> )' --INPUT='%s' --OUTPUT='%s' ",
          algorithm,
          input,
          output
        )
      system(cmd)
      ## 1.3 tif salvestamine tmp kataloogi
      conf$schema <- obj
      conn <- ruut::db_connect(conf = conf)
      grid_layer <-
        sf::read_sf(dsn = conn,
                    layer = c(conf$schema, "bb2_epk2t_grid"))
      DBI::dbDisconnect(conn)
      pk_attributes <- attributes(grid_layer$geom)
      extent <-
        sprintf(
          "%s,%s,%s,%s [EPSG:3301]",
          round((pk_attributes$bbox["xmin"] / 100), digits = 0) * 100,
          ceiling((pk_attributes$bbox["xmax"] / 100)) * 100,
          round((pk_attributes$bbox["ymin"] / 100), digits = 0) * 100,
          ceiling((pk_attributes$bbox["ymax"] / 100)) * 100
        )
      # ruut::qgis_algorithm_search_by_word("Rasterize")
      algorithm <- "gdal:rasterize"
      # ruut::qgis_show_help(algorithm = algorithm)
      input <- vars$tmp_gpkg_file_input_2
      tif_file_name <-
        sprintf("tif_stat_%s", andmed$table[i])
      output <- sprintf("%s/%s.tif", vars$tmp_dir, tif_file_name)
      str <-
        sprintf(
          "{ 'BURN' : 0, 'DATA_TYPE' : 5, 'EXTENT' : '%s', 'EXTRA' : '-a_srs epsg:3301', 'FIELD' : 'value', 'HEIGHT' : 100, 'INIT' : None, 'INPUT' : '%s', 'INVERT' : False, 'NODATA' : -1, 'OPTIONS' : '', 'OUTPUT' : '%s', 'UNITS' : 1, 'WIDTH' : 100 }",
          extent,
          input,
          output
        )
      cmd <-
        ruut::construct_qgis_output_result_to_better_format(str = str, algorithm = algorithm)
      system(cmd) # Salvestame kataloogi
      ## 1.4 ----------- TIF TO POSTGIS -------------
      (
        cmd <- sprintf(
          "export PGPASSWORD=%s && raster2pgsql -s 3301 -d -I -C -M %s %s.%s_tif | psql -U %s -d %s -h %s -p %s",
          conf$password,
          output,
          conf$schema,
          tif_file_name,
          conf$user,
          conf$dbname,
          conf$host,
          conf$port
        )
      )
      system(cmd)
      ## 1.5 Eemaldame üleliigsed veerud ja salvestame tiff alusfaili andmebaasi
      # ruut::qgis_algorithm_search_by_word("Drop ")
      algorithm <- "native:deletecolumn"
      input <- vars$tmp_gpkg_file_input_2
      conf$table <-
        sprintf("grid_stat_%s", andmed$table[i])
      output <-
        ruut::construct_to_gpkg_output_postgres_str(
          conf = conf,
          geometry_field = "geom",
          geometry_type = "LineString",
          srid = 3301,
          checkPrimaryKeyUnicity = FALSE,
          key = "id"
        )
      cmd <-
        sprintf(
          "qgis_process run %s --COLUMN='fid' --COLUMN='left' --COLUMN='top' --COLUMN='right' --COLUMN='bottom' --COLUMN='rahvaarv' --INPUT='%s' --OUTPUT='%s' ",
          algorithm,
          input,
          output
        )
      system(cmd)
      ## 1.6 Kustutame geom veeru.
      conf$schema <- obj
      conn <- ruut::db_connect(conf = conf)
      q <-
        sprintf("ALTER TABLE \"%s\".\"%s\" DROP COLUMN geom;",
                conf$schema,
                conf$table)
      cat(sprintf("\n-----------------\n%s\n\n", q))
      DBI::dbSendQuery(conn, q)
      ## 1.7 Vaakum.
      q <-
        sprintf(
          "VACUUM (FULL, FREEZE, VERBOSE, ANALYZE, INDEX_CLEANUP) \"%s\".\"%s\";",
          conf$schema,
          conf$table
        )
      cat(sprintf("\n-----------------\n%s\n\n", q))
      DBI::dbSendQuery(conn, q)
      ## 1.8 Add foreign key.
      q <-
        sprintf(
          "ALTER TABLE \"%s\".\"%s\" ADD CONSTRAINT %s_fk FOREIGN KEY (id) REFERENCES %s.bb2_epk2t_grid(id);",
          conf$schema,
          conf$table,
          conf$table,
          conf$schema
        )
      cat(sprintf("\n-----------------\n%s\n\n", q))
      DBI::dbSendQuery(conn, q)
      DBI::dbDisconnect(conn)
      ## ----------- 1x1 km rahvaarvu ülekandmine 100x100m ruutudesse ----------
      ## Selleks leiame elamute aadressandmed 1x1km ruutudes, arvutame elanike
      ## keskmise arvu elamutes (s.o elanike arv / majade arv). Edasi leitud
      ## keskmise kaudu kanname elanike arvud 100x100m asuvatesse ruutudesse.
      conf$schema <- obj
      # 2.1 Elamute arv 1x1 km ruutudes
      algorithm <- "native:countpointsinpolygon"
      output <- vars$tmp_gpkg_file_output_1
      conf$table <- "data_p_ma_aadressandmed_EE"
      points <-
        ruut::construct_to_gpkg_output_postgres_str(
          conf = conf,
          geometry_field = "geom",
          geometry_type = NULL,
          srid = 3301,
          checkPrimaryKeyUnicity = FALSE,
          key = "id"
        )
      conf$table <- "bb2_epk2t_grid"
      polygons <-
        ruut::construct_to_gpkg_output_postgres_str(
          conf = conf,
          geometry_field = "geom",
          geometry_type = NULL,
          srid = 3301,
          checkPrimaryKeyUnicity = FALSE,
          key = "id"
        )
      str <- sprintf(
        "{
    'FIELD' : 'numpoints',
    'OUTPUT' : '%s',
    'POINTS' : '%s',
    'POLYGONS' : '%s'}",
        output,
        points,
        polygons
      )
      cmd <-
        ruut::construct_qgis_output_result_to_better_format(str, algorithm)
      system(cmd)
      ## 2.2 tif salvestamine tmp kataloogi
      conn <- ruut::db_connect(conf = conf)
      # 1x1km ruudud
      grid_layer <-
        sf::read_sf(dsn = conn,
                    layer = c(conf$schema, "bb2_epk2t_grid"))
      DBI::dbDisconnect(conn)
      pk_attributes <- attributes(grid_layer$geom)
      extent <-
        sprintf(
          "%s,%s,%s,%s [EPSG:3301]",
          round((pk_attributes$bbox["xmin"] / 100), digits = 0) * 100,
          ceiling((pk_attributes$bbox["xmax"] / 100)) * 100,
          round((pk_attributes$bbox["ymin"] / 100), digits = 0) * 100,
          ceiling((pk_attributes$bbox["ymax"] / 100)) * 100
        )
      # ruut::qgis_algorithm_search_by_word("Rasterize")
      algorithm <- "gdal:rasterize"
      # ruut::qgis_show_help(algorithm = algorithm)
      input <- vars$tmp_gpkg_file_input_1
      tif_file_name <-
        sprintf("tif_stat_%s", "p_ma_aadressandmed_ee")
      output <- sprintf("%s/%s.tif", vars$tmp_dir, tif_file_name)
      str <-
        sprintf(
          "{ 'BURN' : 0, 'DATA_TYPE' : 5, 'EXTENT' : '%s', 'EXTRA' : '-a_srs epsg:3301', 'FIELD' : 'numpoints', 'HEIGHT' : 100, 'INIT' : None, 'INPUT' : '%s', 'INVERT' : False, 'NODATA' : -1, 'OPTIONS' : '', 'OUTPUT' : '%s', 'UNITS' : 1, 'WIDTH' : 100 }",
          extent,
          input,
          output
        )
      cmd <-
        ruut::construct_qgis_output_result_to_better_format(str = str, algorithm = algorithm)
      system(cmd)
      ## 2.3 ----------- TIF TO POSTGIS -------------
      (
        cmd <- sprintf(
          "export PGPASSWORD=%s && raster2pgsql -s 3301 -d -I -C -M %s %s.%s_tif | psql -U %s -d %s -h %s -p %s",
          conf$password,
          output,
          conf$schema,
          tif_file_name,
          conf$user,
          conf$dbname,
          conf$host,
          conf$port
        )
      )
      system(cmd)
      ## 2.4 Elanike tihedus majades. Salvestame TIF-na andmebaasi.
      majade_arv <- "tif_stat_p_ma_aadressandmed_ee_tif"
      rahvaarv <- sprintf("tif_stat_%s_tif", andmed$table[i])
      tihedus <- sprintf("tif_stat_density_%s_tif", andmed$table[i])
      # F-n tiheduse leidmiseks mapalgebra f-ni jaoks
      conn <- ruut::db_connect(conf = conf)
      q <- sprintf(
        "create or replace function
  rahvaarv_tihedus_epk2t_fn(pixel float[][][], pos integer[][], variadic userargs text[])
  returns float
  language plpgsql
  immutable  -- careful: this function is immutable, yours may not be
as $$
  declare
    pixval1 float;
    pixval2 float;
    tihedus float;
  begin
    pixval1 := pixel[1][1][1];  -- pixel indices: [raster #][xdistance][ydistance]
    pixval2 := pixel[2][1][1];  -- pixel indices: [raster #][xdistance][ydistance]
    if pixval1 > 0  then
      tihedus := (pixval2 / pixval1);
    else
      tihedus := 0;
    end if;
    return tihedus;
  end;
$$;"
      )
      cat(sprintf("\n-----------------\n%s\n\n", q))
      DBI::dbSendQuery(conn, q)
      # Salvestame andmebaasi
      q <- sprintf("
DROP TABLE IF EXISTS  \"%s\".\"%s\";
",
                   conf$schema,
                   tihedus)
      cat(sprintf("\n-----------------\n%s\n\n", q))
      DBI::dbSendQuery(conn, q)
      q <-
        sprintf(
          "SELECT ST_MapAlgebra(t2.rast, 1, t1.rast, 1, 'rahvaarv_tihedus_epk2t_fn(double precision[], int[], text[])'::regprocedure) AS rast INTO \"%s\".\"%s\" FROM \"%s\".\"%s\" t1,  \"%s\".\"%s\" t2;",
          conf$schema,
          tihedus,
          conf$schema,
          rahvaarv,
          conf$schema,
          majade_arv
        )
      cat(sprintf("\n-----------------\n%s\n\n", q))
      DBI::dbGetQuery(conn, q)
      DBI::dbDisconnect(conn)
      ## 2.2 Rahvaarv 100x100m ruutudes.
      tihedus <- sprintf("tif_stat_density_%s_tif", andmed$table[i])
      majade_arv <- "tif_stat_p_ma_aadressandmed_ee_tif"
      rahvaarv <- uus_tabel
      conn <- ruut::db_connect(conf = conf)
      q <- sprintf("
DROP TABLE IF EXISTS  \"%s\".\"%s\";
",
                   conf$schema,
                   rahvaarv)
      cat(sprintf("\n-----------------\n%s\n\n", q))
      DBI::dbSendQuery(conn, q)
      q <- sprintf(
        "SELECT
    ST_MapAlgebra(
        t1.rast, 1,
        t2.rast, 1,
        '([rast2.val] * [rast1.val])'
    ) AS rast INTO \"%s\".tif_stat_rel_1x1km_epk02t_tif
FROM \"%s\".tif_p_ma_aadressandmed_ee_tif t1
cross JOIN \"%s\".tif_stat_density_rel_1x1km_tif t2;",
        conf$schema,
        conf$schema,
        conf$schema
      )
      cat(sprintf("\n-----------------\n%s\n\n", q))
      DBI::dbGetQuery(conn, q)
      DBI::dbDisconnect(conn)
    }
  }
  ## Layers list
  conf$schema <- obj
  ruut::db_schema_tablenames(conf = conf)
  sf::st_layers(dsn = vars$tmp_gpkg_file)
}