Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- # make dead-ends ratio map
- library('ks')
- library('raster')
- # connect to database
- library("RPostgreSQL")
- driver <- dbDriver("PostgreSQL")
- connection <- dbConnect(driver, dbname="", password="", user="", host="localhost")
- cities = read.csv('cities.csv',stringsAsFactors=FALSE)
- # iterate over cities
- for( i in 1:length(cities$name)){
- city = cities[i,'name']
- epsg = cities[i,'epsg']
- print(paste('working on',city))
- d = dbGetQuery(connection,paste(
- "SELECT",
- "km,flags,",
- "car_comp AS cc,",
- "bike_comp AS bc,",
- "foot_comp AS fc,",
- "car_directness AS cd,",
- "bike_directness AS bd,",
- "foot_directness AS fd,",
- "ST_X(ST_Transform(ST_Centroid(geom_way),",epsg,")) AS x,", # meters
- "ST_Y(ST_Transform(ST_Centroid(geom_way),",epsg,")) AS y", # meters
- "FROM",city
- ))
- # calculate grid resolution
- cell_size = 250 # square cell size in METERS
- bandwidth = 2000^2 # meters^2
- xres = round((max(d$x)-min(d$x)) / cell_size)
- yres = round((max(d$y)-min(d$y)) / cell_size)
- # subset for each mode
- car = d[!is.na(d$cc),c('x','y','km','cc','cd')]
- #bike = d[!is.na(d$bc),c('x','y','km','bc','bd')]
- #foot = d[!is.na(d$fc),c('x','y','km','fc','fd')]
- # calculate the KDE
- total_kde = kde(
- x=matrix(c(car$x,car$y),length(car$x),2),
- w=car$km / sum(car$km) * length(car$km),
- H=matrix(c(bandwidth,0,0,bandwidth),2),
- xmin=c(min(d$x),min(d$y)),
- xmax=c(max(d$x),max(d$y)),
- gridsize=c(xres,yres),
- approx.cont=FALSE,
- verbose=TRUE
- )
- # rasterize
- total = raster(total_kde)
- projection(total) = CRS(paste0('+init=epsg:',epsg))
- # scale to the length of road
- total = (total / sum(total_kde$estimate)) * sum(car$km)
- # deadends
- deadcar = car[car$cc==-1,]
- dead_kde = kde(
- x=matrix(c(deadcar$x,deadcar$y),length(deadcar$x),2),
- w=deadcar$km / sum(deadcar$km) * length(deadcar$km),
- H=total_kde$H, # reuse H from the total version
- xmin=c(min(d$x),min(d$y)),
- xmax=c(max(d$x),max(d$y)),
- gridsize=c(xres,yres),
- verbose=TRUE
- )
- # rasterize
- dead = raster(dead_kde)
- projection(dead) = CRS(paste0('+init=epsg:',epsg))
- # scale to the length of dead-end road
- dead = (dead / sum(dead_kde$estimate)) * sum(deadcar$km)
- # Percent Dead-ends
- percent = dead / total
- percent[is.na(percent)] = 0
- percent[percent>1] = 1
- banded_raster = addLayer(total,dead,percent)
- names(banded_raster) = c('car_total','car_dead','car_percent')
- # output the file as geotiff
- writeRaster(
- banded_raster,
- filename=paste0('data/ratio-rasters/',city,'_car.tiff'),
- format='GTiff',
- overwrite=TRUE
- )
- }
- dbDisconnect(connection)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement