#' 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( )' --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) }