https://docs.google.com/document/d/13Q6PfnHiiwitlz87VkDdADERi3MXKRRs17igEVppCLk/edit#

R setup

# clear R memory
rm(list = ls())

# db_connect.r
source("P:/db_connect.r")

# import libraries
library(tidyverse)
library(sf)
library(lwgeom)
library(geojsonio)

options(scipen = 1000000)
options(digits = 6)

#sessionInfo()

dat

```{sql connection=statop, output.var = "einw"}

select BFS_NR, GEMEINDE, count(*) as count from DATASHOP.EINW_GEOKOD_HIST where stichtag = to_date('30.09.2019', 'dd.mm.yyyy') and zivilrechtlich = 1 group by BFS_NR, GEMEINDE

# create function topojson_write_stat

```r

geodat_path ="C:/gitrepos/Ressources_Maps/2019/GemeindeGrosseSeeOhneExklave_gen_epsg2056_F_KTZH_2019.json"
geodat_bfsnr_char = "GDE_ID";
dat_df = einw;
dat_bfsnr_char = "BFS_NR";
topojson_path = "C:/gitrepos/Ressources_Maps/2019/GemeindeGrosseSeeOhneExklave_gen_epsg2056_F_KTZH_2019.topojson";
pixel_width_int = 480

# Funktion topojson_write_stat
topojson_write_stat <- function(geodat_path, geodat_bfsnr_char, dat_df, dat_bfsnr_char, topojson_path, pixel_width_int) {

  # import geodat
  geodat <- sf::st_read(geodat_path, stringsAsFactors = FALSE, crs=2056)

  # geometry: edit scaling
  #print(st_bbox(geodat))
  bbox <- as.integer(st_bbox(geodat))
  scaling <- pixel_width_int/(bbox[3]-bbox[1])
  # 1) scaling 
  geodat_sfc <- st_geometry(geodat)*scaling 
  print(st_bbox(geodat_sfc))
  # 2) subtract bbox 
  geodat_sfc <- st_geometry(geodat)*scaling-c(bbox[1]*scaling, bbox[2]*scaling)
  print(st_bbox(geodat_sfc))
  bbox2 <- as.integer(st_bbox(geodat_sfc))
  # 3) snap to grid
  geodat_sfc2 <- st_snap_to_grid(geodat_sfc, 1) # coordinates don't have decimal places anymore
  print(st_bbox(geodat_sfc2))
  # View(geodat_sf3)
  # 4) rescale y axis: y1 = -1*(y-yMax) (drehen der y koordinate)
  # rescale y axis part 1:  y-yMax
  geodat_sfc3 <- geodat_sfc2 - c(0,bbox2[4])
  print(st_bbox(geodat_sfc3))
  # rescale y axis part 2: *-1
  # z <- geodat_sf1[1,]
  ## function
  geodat_sfc4 <- geodat_sfc3 * matrix(data = c(1,-1), ncol = 2)
  st_coordinates(geodat_sfc4)
  print(st_bbox(geodat_sfc4))
  #View(geodat_sf2)
  #View(st_coordinates(geodat_sf1)); View(st_coordinates(geodat_sf2)) 
  geodat_sf4 <- st_sf(st_geometry(geodat_sfc4), geodat %>% st_drop_geometry())
  View(geodat_sf4)
  print(st_bbox(geodat_sf4))

  # attribute: join sql dat  
  # https://stackoverflow.com/questions/54823846/dplyr-left-join-does-not-work-with-a-character-objects-as-the-lhs-variable
  join_cols <- c(dat_bfsnr_char)
  names(join_cols) <- geodat_bfsnr_char
  geodat_attr <- geodat_sf4 %>% st_drop_geometry() %>% left_join(dat_df,by = join_cols)
  #print("geodat_attr is created")


  # check
  if(
    nrow(geodat) != nrow(geodat_attr)
  )
  {stop (paste("nrow(geodat) should be == nrow(dat)"))} 
  print("ok: nrow(geodat) == nrow(dat)")

  # sf object
  geodat_sf = st_sf(geodat_attr,st_geometry(geodat_sf4)) 
  print("geodat_sf is created")

  # write topojson
  geojsonio::topojson_write(
  geodat_sf, 
  lat = NULL, 
  lon = NULL,
  geodatetry = "polygon",
  group = NULL,
  file = topojson_path,
  overwrite = TRUE,
  precision =  NULL,
  convert_wgs84 = FALSE, 
  crs = 2056,
  object_name = "test", 
  quantization = 0)
}

Run function topojson_write_stat

# L:/STAT/08_DS/03_GIS/Geodaten/2019/WahlkreiseAuszaehlkreiseGrosseSeeOhneExklaveSitze_gen_F_KTZH_2019.shp => hat funktioniert

# @geodat_path: geodata input path
# @geodat_bfsnr_char: geodat-id by which geodat should be joined to dat
# @dat_df: data that should be joined to geodat. 
# => dat must have an id that matches geodat. See: dat_bfsnr_char
# => dat must be aggragated to match geodat. E.g. when geodat = community then dat has to have only one row for each community. 
# @dat_bfsnr_char:  dat-id by which dat should be joined to geodat
# @topojson_path: topojson output path

topojson_write_stat(
  geodat_path = "L:/STAT/08_DS/03_GIS/Geodaten/2019/GemeindeGrosseSeeOhneExklave_gen_F_KTZH_2019.shp",
                    geodat_bfsnr_char = "GDE_ID", 
                    dat_df = einw, 
                    dat_bfsnr_char = "BFS_NR", 
                    topojson_path = "C:/gitrepos/Ressources_Maps/2019/GemeindeGrosseSeeOhneExklave_gen_epsg2056_F_KTZH_2019.topojson", 
                    pixel_width_int = 480)


statistikZH/Ressources_Maps documentation built on Dec. 23, 2021, 5:28 a.m.