210127_artikkel/maatriksid.R

## Maatriksid
# 0 - intersect
# 1 - contain
# 2 - disjoint
# 3 - equal
# 4 - touch
# 5 - overlap
# 6 - are within
# 7 - cross


source("functions/delete_existing_variables.R")
library(dplyr)
library(sf)

gpkg_home <- "/data/gpkg/artiklid/artikkel_210127_valga_matsalu_lahemaa"
obj <- "marja"
dsn <- sprintf("%s/%s.gpkg", gpkg_home, obj)
conf <- ruut::get_config()
conf$gpkg_home <- gpkg_home
conf$gpkg_file <- obj

ruut::gpkg_andmebaasi_kihtide_nimekiri(obj = obj, gpkg_home = gpkg_home)



## Kõik geomeetrilised objektid GPKG andmebaasist
## Loeme andmebaasist piiri ja piirikasti.
# Layers list
gpkg_info <- sf::st_layers(dsn = dsn)
layer_names <- gpkg_info$name
for (layer_name in layer_names) {
  cat(sprintf("\n%s", layer_name))
  assign(layer_name, sf::read_sf(dsn = dsn, layer = layer_name))
}




## 0-maatriks
# bb_epk02t_grid <- sf::read_sf(dsn = dsn, layer = "bb_epk02t_grid")
rows <- sort(unique(bb_epk02t_grid$bottom), decreasing = T)
cols <- sort(unique(bb_epk02t_grid$left))
length(rows) * length(cols)

## Nullmaatriks
(m.0 <- matrix(0L, nrow = length(rows), ncol = length(cols), byrow = F)) # 0-maatriks
## Ruudu id väärtustega maatriks
(m.id <- matrix(bb_epk02t_grid$id, nrow = length(rows), ncol = length(cols), byrow = F))
## Ruudu vasakpoolse koordinaadi väärtustega maatriks
(m.left <- matrix(bb_epk02t_grid$left, nrow = length(rows), ncol = length(cols), byrow = F))
## Ruudu parempoolse koordinaadi väärtustega maatriks
(m.right <- matrix(bb_epk02t_grid$right, nrow = length(rows), ncol = length(cols), byrow = F))
## Ruudu ülemise koordinaadi väärtustega maatriks
(m.top <- matrix(bb_epk02t_grid$top, nrow = length(rows), ncol = length(cols), byrow = F))
## Ruudu alumise koordinaadi väärtustega maatriks
(m.bottom <- matrix(bb_epk02t_grid$bottom, nrow = length(rows), ncol = length(cols), byrow = F))

## ----------------------- TRUE/FALSE matrix ----------------------
## Kas alusruudustik sisaldab meie valitud kihti?

x <- layer_names[11]
x <- "piir"
x <- "bb_e_401_hoone_ka_10"
x <- "bb_muutee"
x <- "bb_aadressandmed"


for (x in layer_names) {
  ruumiline_obj <- get(x)
  y <- unlist(sf::st_intersects(ruumiline_obj, bb_epk02t_grid, sparse = TRUE))
  z <- rep(0, length(bb_epk02t_grid$id))
  z[y] <- 1
  assign(sprintf("m.%s", x), matrix(z, nrow = length(rows), ncol = length(cols), byrow = F))
  get(sprintf("m.%s", x))

  ruudustik <- bb_epk02t_grid
  ruudustik$value <- z

  ## plot to file
  png(filename = sprintf("tmp/img/matrix_true_false/%s.png", x))
  sf::st_geometry(ruudustik) %>% plot(main = x, sub = "True/False")
  sf::st_geometry(ruumiline_obj) %>% plot(add = T, border = "dark red", lwd = 1, col = "blue", pch = 16)
  text(x = (ruudustik$left + ruudustik$right) / 2, y = (ruudustik$bottom + ruudustik$top) / 2, labels = as.character(ruudustik$value))
  sf::st_geometry(ruudustik) %>% plot(add = T)
  dev.off()


  ## Punktide arv ruudus
  ## Kontlrollime kas geomeetriline objekt on punkt.
  (is_point <- any(grepl("point", tolower(attributes(ruumiline_obj$geom)$class), fixed = TRUE)))
  if (is_point) {
    # ruut::qgis_algorithm_search_by_word("Count")
    algorithm <- "native:countpointsinpolygon"
    conf$gpkg_table <- sprintf("%s_numpoints", x)
    output <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = F)
    conf$gpkg_table <- sprintf("%s", x)
    points <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = T)
    conf$gpkg_table <- sprintf("%s", "bb_epk02t_grid")
    polygons <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = T)
    str <- sprintf("{ 'CLASSFIELD' : '', 'FIELD' : 'value', 'OUTPUT' : '%s', 'POINTS' : '%s', 'POLYGONS' : '%s', 'WEIGHT' : '' }", output, points, polygons)
    cmd <- ruut::construct_qgis_output_result_to_better_format(str = str, algorithm = algorithm)
    system(cmd)
  }



  ## Kontlrollime kas geomeetriline objekt on polügoon.
  ## Arvutame sel juhul pindala.
  is_polygon <- any(grepl("polygon", tolower(attributes(ruumiline_obj$geom)$class), fixed = TRUE))
  is_line <- any(grepl("line", tolower(attributes(ruumiline_obj$geom)$class), fixed = TRUE))

  if (is_polygon || is_line) {
    ## Lõikame objekti ruudistikga tükkideks. Ühte ruutu võib jääda mitu tükki.
    ## Peame need pärast ühendama
    # ruut::qgis_algorithm_search_by_word("Intersection")
    algorithm <- "native:intersection"
    conf$gpkg_table <- sprintf("%s", x)
    input <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = T)
    conf$gpkg_table <- sprintf("%s", "test_layer")
    output <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = F)
    conf$gpkg_table <- sprintf("%s", "bb_epk02t_grid")
    overlay <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = T)
    str <- sprintf("{ 'INPUT' : '%s', 'INPUT_FIELDS' : ['fid'], 'OUTPUT' : '%s', 'OVERLAY' : '%s', 'OVERLAY_FIELDS' : [], 'OVERLAY_FIELDS_PREFIX' : '' }", input, output, overlay)
    cmd <- ruut::construct_qgis_output_result_to_better_format(str = str, algorithm = algorithm)
    system(cmd)
    ## Ühendame ruudus olevad pinnad
    # ruut::qgis_algorithm_search_by_word("Dissolve")
    algorithm <- "native:dissolve"
    conf$gpkg_table <- sprintf("%s", "test_layer")
    input <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = T)
    conf$gpkg_table <- sprintf("%s", "test_layer_2")
    output <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = F)
    str <- sprintf("{ 'FIELD' : ['%s'], 'INPUT' : '%s', 'OUTPUT' : '%s' }", "id", input, output)
    cmd <- ruut::construct_qgis_output_result_to_better_format(str = str, algorithm = algorithm)
    system(cmd)
    # processing.runalg('qgis:addfieldtoattributestable', input_layer, field_name, field_type, field_length, field_precision, output_layer)
      # processing.runalg('qgis:advancedpythonfieldcalculator', input_layer, field_name, field_type, field_length, field_precision, global, formula, output_layer)
    if (is_polygon) {
      ## Lisame pindade pindalad ja perimeetri
      # ruut::qgis_algorithm_search_by_word("attributes")
      algorithm <- "qgis:exportaddgeometrycolumns"
      conf$gpkg_table <- sprintf("%s", "test_layer_2")
      input <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = T)
      conf$gpkg_table <- sprintf("%s", "test_layer")
      output <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = F)
      str <- sprintf("{ 'CALC_METHOD' : 0, 'INPUT' : '%s', 'OUTPUT' : '%s' }", input, output)
      cmd <- ruut::construct_qgis_output_result_to_better_format(str = str, algorithm = algorithm)
      system(cmd)
      ## Add field
      ## Arvutame pindala protsendi ruudust
      # ruut::qgis_algorithm_search_by_word("pythonfield")
      algorithm <- "qgis:advancedpythonfieldcalculator"
      conf$gpkg_table <- sprintf("%s", "test_layer")
      input <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = T)
      conf$gpkg_table <- sprintf("%s_area", x)
      output <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = F)
      str <- sprintf("{ 'FIELD_NAME' : 'area_pc', 'FIELD_TYPE' : 1, 'FIELD_LENGTH' : 4, 'FIELD_PRECISION' : 2, 'GLOBAL' : '', 'FORMULA' : 'value = round($geom.area()/10001,4)', 'INPUT' : '%s', 'OUTPUT' : '%s' }", input, output)
      cmd <- ruut::construct_qgis_output_result_to_better_format(str = str, algorithm = algorithm)
      system(cmd)
    } else {
      ## Lisame joonte pikkused
      # ruut::qgis_algorithm_search_by_word("attributes")
      algorithm <- "qgis:exportaddgeometrycolumns"
      conf$gpkg_table <- sprintf("%s", "test_layer_2")
      input <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = T)
      conf$gpkg_table <- sprintf("%s_length", x)
      output <- ruut::construct_to_gpkg_output_file_str(conf = conf, geometry_field = "geom", is_input_str = F)
      str <- sprintf("{ 'CALC_METHOD' : 0, 'INPUT' : '%s', 'OUTPUT' : '%s' }", input, output)
      cmd <- ruut::construct_qgis_output_result_to_better_format(str = str, algorithm = algorithm)
      system(cmd)
    }






  }
}















## ----------------- Confusion matrix -------------------
# \url{https://en.wikipedia.org/wiki/Confusion_matrix}


## ----------------- export matrix -------------------
mat <- matrix(1:10, ncol = 2)
rownames(mat) <- letters[1:5]
colnames(mat) <- LETTERS[1:2]

mat
write.table(mat, file = "test.txt") # keeps the rownames
read.table("test.txt", header = TRUE, row.names = 1) # says first column are rownames
unlink("test.txt")
write.table(mat, file = "test2.txt", row.names = FALSE) # drops the rownames
read.table("test.txt", header = TRUE)
unlink("test2.txt")