In MariekeDirk/uhi_max: Maximum Urban Heat island
knitr::opts_chunk$set(echo = TRUE)
library(uhimax)
library(raster)
library(rgdal)
library(mapview)
library(data.table)
library(tidyverse)
library(lubridate)
library(dplyr)
library(caret)
Initialization
The urban heat island equations are based on Theeuwes (2017).
wd<-"D:/uhimax/"
grn<-raster("D:/groenkaart_rivm/rivm_20170415_g_groenkaart_10m.tif")
# grd.svf<-raster(paste0(wd,"Grids_veg_svf/svf_utrecht_1m.grd"))
# grd.svf<-blur_svf_fveg(grd.svf)
# grd.svf<-projectRaster(grd.svf,crs=crs(grn))
#
# grn.fveg<-crop(grn,extent(grd.svf))
# grn.fveg<-blur_svf_fveg(grn.fveg)/100
# grd.fveg<-mask(grn.fveg,grd.svf)
#
# st_fveg_svf<-stack(grd.svf,grd.fveg)
# names(st_fveg_svf)<-c("svf","fveg")
# writeRaster(st_fveg_svf,filename=paste0(wd,"Grids_veg_svf/svf_grn_blur500m.grd"),
# overwrite=TRUE)
#
# mapview(grd.fveg)
Wunderground network data in Utrecht
Within the city of Utrecht we have carefully selected three stations from the Wunderground network:
- IUTRECHT23 (Used by Theeuwes (2017))
- IUTRECHT196 (not in sensitivity area)
- IUTRECHT376
- IUTRECHT299
These stations have been measuring the city temperature for the past couple of years. The data can be downloaded using download_time_seq, note that the temperature is in Farenheid and not degrees Celcius. The data has been filtered and interpolated to 10minute timestamps. Unrealistic high and low temperatures are set to NA. In case more than 8 identical measurements in a row occur the value is also set to NA.
svf_grn<-readRDS("D:/uhimax/Wunderground/wunderground.rds")
df.svf_grn<-data.frame(svf_grn)
stations<-df.svf_grn$Station.ID
city_T<-paste0(wd,"Wunderground/Filtered/",df.svf_grn$Station.ID,"_filtered.rds")
svf<-df.svf_grn$svf
fveg<-df.svf_grn$fveg
sp.IUTRECHT<-svf_grn
# coordinates(sp.IUTRECHT)<-~Lon+Lat
# crs(sp.IUTRECHT)<-crs("+init=epsg:28992")
From 10min data to Daily UHImax parameters
From the autmaomatic weather stations (AWS) in the rural area the following data is required:
- 10 minute global radiation
- 10 minute wind speed at 10m
- 10 minute temperature at 1.5m
- 10 minute pressure derived to MSL
- 1 hour relative humidity
- 1 hour precipitation
- 1 hour mean wind
The 10 minute data from global radiation, wind speed, temperature and pressure is used to calculate the UHImax. The hourly relative humidity, precipitation and mean wind are used to filter out synoptic situations with frontal system and fog. Under these different conditions there is no clear relation between the AWS and city temperature.
Cabauw_data<-prepare_Cabauw(ten_min = paste0(wd,"Cabauw_meteo/cabauw10min.csv"),
hourly_data = paste0(wd,"Cabauw_meteo/cabauw_hourly.csv"))
Calculate UHI max parameters
Setting all parameters to required to calculate the UHI with the formula of Theeuwes (2017). For this figure no subset of meteorological conditions was made.
# uhi_Utrecht<-uhimax_params(STN="UTRECHT23",
# svf=IUtrecht23_svf,
# fveg=IUtrecht23_fveg,
# city_T=paste0(wd,"Wunderground/Filtered/IUTRECHT23_filtered.rds"),
# rural_T=Cabauw_data$Cabauw_10min,
# rural_meteo=Cabauw_data$Cabauw_Theeuwes)
# ggplot(uhi_Utrecht,aes(Tref,Tcity,colour=S_new)) + geom_point() + geom_abline() + scale_color_gradientn(colours = heat.colors(20))
#Similarly the same function can be applied to other stations within Utrecht:
Wunderground_Utrecht<-mapply(uhimax_params, STN = stations, svf = svf, fveg = fveg, city_T = city_T,
MoreArgs = list(
rural_T = Cabauw_data$Cabauw_10min,
rural_meteo = Cabauw_data$Cabauw_Theeuwes
#,wd = paste0(wd,"UHImax/"), write.file = TRUE
),SIMPLIFY = FALSE)
Wunderground_Utrecht<-do.call("rbind",Wunderground_Utrecht)
ggplot(Wunderground_Utrecht,aes(Tref,Tcity,colour=S_new)) + geom_point() + geom_abline() + scale_color_gradientn(colours = heat.colors(20))
Subset based on weather conditions
uhimax_files<-list.files(paste0(wd,"UHImax/"),
full.names = TRUE,pattern=".txt")
uhimax_files<-lapply(uhimax_files,fread)
uhimax_Utrecht<-do.call("rbind",uhimax_files)
uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$stn %in% c("UTRECHT23","IUTRECHT376","IUTRECHT299")),]
# uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$stn %in% c("IUTRECHT299")),]
# uhimax_Utrecht$stn[which(uhimax_Utrecht$stn=="UTRECHT23")]<-"IUTRECHT23"
# uhimax_Utrecht<-fread(paste0(wd,"UHImax/IUTRECHT23_UHIparams.txt"))
uhimax_Utrecht<-merge(Cabauw_data$Cabauw_sub,uhimax_Utrecht,by=c("start","stop"))
uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Rain==TRUE),]
uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Wind==TRUE),]
uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$rh==TRUE),]
# uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Select==TRUE),]
uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Tref>17),]
# uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$stn %in% c("IUTRECHT23","IUTRECHT376")),]
wur.R2<-caret::R2((2-uhimax_Utrecht$svf-uhimax_Utrecht$fveg)*uhimax_Utrecht$meteo,uhimax_Utrecht$UHImeasured)
wur.RMSE<-caret::RMSE((2-uhimax_Utrecht$svf-uhimax_Utrecht$fveg)*uhimax_Utrecht$meteo,uhimax_Utrecht$UHImeasured)
Correlations with the three meteo params
The maximum urban heat island is in the cases we selected related to:
$$ UHI_{max} = (2-SVF-F_{veg}) \sqrt[4]{\frac{S*DTR^3}{U}} $$
The station IUTRECHT23 has the best correlation coefficient with 0.6, while also including the other observations in the center of Utrecht reduces R2 to 0.4.
uhi_u<-ggplot(uhimax_Utrecht,aes(U,UHImeasured,colour=DTR)) +geom_point() + scale_color_gradientn(colours = terrain.colors(20))
uhi_DTR<-ggplot(uhimax_Utrecht,aes(DTR,UHImeasured,colour=S_new)) + geom_point()+scale_color_gradientn(colours = heat.colors(20))
uhi_S<-ggplot(uhimax_Utrecht,aes(S_new,UHImeasured,colour=U)) +geom_point() + geom_point()+scale_color_gradientn(colours = topo.colors(20))
ggsave(uhi_u,filename = paste0(wd,"UHImax/fig/uhi_u.png"))
ggsave(uhi_DTR,filename = paste0(wd,"UHImax/fig/uhi_DTR.png"))
ggsave(uhi_S,filename = paste0(wd,"UHImax/fig/uhi_S.png"))
uhi_u
uhi_DTR
uhi_S
caret::R2((2-uhimax_Utrecht$svf-uhimax_Utrecht$fveg)*uhimax_Utrecht$meteo,uhimax_Utrecht$UHImeasured)
caret::RMSE((2-uhimax_Utrecht$svf-uhimax_Utrecht$fveg)*uhimax_Utrecht$meteo,uhimax_Utrecht$UHImeasured)
R2.exp <- expression(paste(" ",R^2 ,"= 0.38"))
uhimax_corr<-ggplot(uhimax_Utrecht,aes(UHImeasured,(2-svf-fveg)*meteo,colour=factor(stn))) + geom_point(size=3) + geom_abline() +
xlim(0,10) + ylim(0,10) + coord_fixed() +
xlab("UHImax measured") + ylab("UHImax calculated") +
scale_color_manual(values=c('orange','yellow','green' )) +
guides(color=guide_legend("Station ID"))+
annotate("text",x=1,y=10,label=R2.exp) +
annotate("text",x=1,y=9.4,label=paste0("RMSE= ",round(wur.RMSE,2)))
uhimax_corr
ggsave(uhimax_corr,
filename="C:/Users/marie/OneDrive/Afbeeldingen/wunderground_utrecht.png")
Using the formula of Theeuwes for Utrecht
The formula of Theeuwes(2017) was derived for:
- 0 < fveg < 0.4
- 0.2 < SVF < 0.9
Outside of this range the formula has not been tested.
svf.r<-raster::resample(grd.svf,grd.fveg,method="bilinear")
values(svf.r)[values(svf.r)<0.2 | values(svf.r)>0.9] = NA
values(grd.fveg)[values(grd.fveg)<0 | values(grd.fveg)>0.4] = NA
city_center<-extent(c(5.1,5.13,52.08,52.098))
city_center<-raster(city_center)
values(city_center)<-1
crs(city_center)<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
city_center<-projectRaster(city_center,crs=crs(grd.fveg))
st<-(2-svf.r-grd.fveg)*mean(uhimax_Utrecht$meteo)
st<-crop(st,city_center)
mapview(st)
MariekeDirk/uhi_max documentation built on May 23, 2019, 1:44 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE)
library(uhimax) library(raster) library(rgdal) library(mapview) library(data.table) library(tidyverse) library(lubridate) library(dplyr) library(caret)
Initialization
The urban heat island equations are based on Theeuwes (2017).
wd<-"D:/uhimax/" grn<-raster("D:/groenkaart_rivm/rivm_20170415_g_groenkaart_10m.tif") # grd.svf<-raster(paste0(wd,"Grids_veg_svf/svf_utrecht_1m.grd")) # grd.svf<-blur_svf_fveg(grd.svf) # grd.svf<-projectRaster(grd.svf,crs=crs(grn)) # # grn.fveg<-crop(grn,extent(grd.svf)) # grn.fveg<-blur_svf_fveg(grn.fveg)/100 # grd.fveg<-mask(grn.fveg,grd.svf) # # st_fveg_svf<-stack(grd.svf,grd.fveg) # names(st_fveg_svf)<-c("svf","fveg") # writeRaster(st_fveg_svf,filename=paste0(wd,"Grids_veg_svf/svf_grn_blur500m.grd"), # overwrite=TRUE) # # mapview(grd.fveg)
Wunderground network data in Utrecht
Within the city of Utrecht we have carefully selected three stations from the Wunderground network:
- IUTRECHT23 (Used by Theeuwes (2017))
- IUTRECHT196 (not in sensitivity area)
- IUTRECHT376
- IUTRECHT299
These stations have been measuring the city temperature for the past couple of years. The data can be downloaded using download_time_seq, note that the temperature is in Farenheid and not degrees Celcius. The data has been filtered and interpolated to 10minute timestamps. Unrealistic high and low temperatures are set to NA. In case more than 8 identical measurements in a row occur the value is also set to NA.
svf_grn<-readRDS("D:/uhimax/Wunderground/wunderground.rds") df.svf_grn<-data.frame(svf_grn) stations<-df.svf_grn$Station.ID city_T<-paste0(wd,"Wunderground/Filtered/",df.svf_grn$Station.ID,"_filtered.rds") svf<-df.svf_grn$svf fveg<-df.svf_grn$fveg sp.IUTRECHT<-svf_grn # coordinates(sp.IUTRECHT)<-~Lon+Lat # crs(sp.IUTRECHT)<-crs("+init=epsg:28992")
From 10min data to Daily UHImax parameters
From the autmaomatic weather stations (AWS) in the rural area the following data is required:
- 10 minute global radiation
- 10 minute wind speed at 10m
- 10 minute temperature at 1.5m
- 10 minute pressure derived to MSL
- 1 hour relative humidity
- 1 hour precipitation
- 1 hour mean wind
The 10 minute data from global radiation, wind speed, temperature and pressure is used to calculate the UHImax. The hourly relative humidity, precipitation and mean wind are used to filter out synoptic situations with frontal system and fog. Under these different conditions there is no clear relation between the AWS and city temperature.
Cabauw_data<-prepare_Cabauw(ten_min = paste0(wd,"Cabauw_meteo/cabauw10min.csv"), hourly_data = paste0(wd,"Cabauw_meteo/cabauw_hourly.csv"))
Calculate UHI max parameters
Setting all parameters to required to calculate the UHI with the formula of Theeuwes (2017). For this figure no subset of meteorological conditions was made.
# uhi_Utrecht<-uhimax_params(STN="UTRECHT23", # svf=IUtrecht23_svf, # fveg=IUtrecht23_fveg, # city_T=paste0(wd,"Wunderground/Filtered/IUTRECHT23_filtered.rds"), # rural_T=Cabauw_data$Cabauw_10min, # rural_meteo=Cabauw_data$Cabauw_Theeuwes) # ggplot(uhi_Utrecht,aes(Tref,Tcity,colour=S_new)) + geom_point() + geom_abline() + scale_color_gradientn(colours = heat.colors(20)) #Similarly the same function can be applied to other stations within Utrecht: Wunderground_Utrecht<-mapply(uhimax_params, STN = stations, svf = svf, fveg = fveg, city_T = city_T, MoreArgs = list( rural_T = Cabauw_data$Cabauw_10min, rural_meteo = Cabauw_data$Cabauw_Theeuwes #,wd = paste0(wd,"UHImax/"), write.file = TRUE ),SIMPLIFY = FALSE) Wunderground_Utrecht<-do.call("rbind",Wunderground_Utrecht) ggplot(Wunderground_Utrecht,aes(Tref,Tcity,colour=S_new)) + geom_point() + geom_abline() + scale_color_gradientn(colours = heat.colors(20))
Subset based on weather conditions
uhimax_files<-list.files(paste0(wd,"UHImax/"), full.names = TRUE,pattern=".txt") uhimax_files<-lapply(uhimax_files,fread) uhimax_Utrecht<-do.call("rbind",uhimax_files) uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$stn %in% c("UTRECHT23","IUTRECHT376","IUTRECHT299")),] # uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$stn %in% c("IUTRECHT299")),] # uhimax_Utrecht$stn[which(uhimax_Utrecht$stn=="UTRECHT23")]<-"IUTRECHT23" # uhimax_Utrecht<-fread(paste0(wd,"UHImax/IUTRECHT23_UHIparams.txt")) uhimax_Utrecht<-merge(Cabauw_data$Cabauw_sub,uhimax_Utrecht,by=c("start","stop")) uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Rain==TRUE),] uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Wind==TRUE),] uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$rh==TRUE),] # uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Select==TRUE),] uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Tref>17),] # uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$stn %in% c("IUTRECHT23","IUTRECHT376")),] wur.R2<-caret::R2((2-uhimax_Utrecht$svf-uhimax_Utrecht$fveg)*uhimax_Utrecht$meteo,uhimax_Utrecht$UHImeasured) wur.RMSE<-caret::RMSE((2-uhimax_Utrecht$svf-uhimax_Utrecht$fveg)*uhimax_Utrecht$meteo,uhimax_Utrecht$UHImeasured)
Correlations with the three meteo params
The maximum urban heat island is in the cases we selected related to:
$$ UHI_{max} = (2-SVF-F_{veg}) \sqrt[4]{\frac{S*DTR^3}{U}} $$
The station IUTRECHT23 has the best correlation coefficient with 0.6, while also including the other observations in the center of Utrecht reduces R2 to 0.4.
uhi_u<-ggplot(uhimax_Utrecht,aes(U,UHImeasured,colour=DTR)) +geom_point() + scale_color_gradientn(colours = terrain.colors(20)) uhi_DTR<-ggplot(uhimax_Utrecht,aes(DTR,UHImeasured,colour=S_new)) + geom_point()+scale_color_gradientn(colours = heat.colors(20)) uhi_S<-ggplot(uhimax_Utrecht,aes(S_new,UHImeasured,colour=U)) +geom_point() + geom_point()+scale_color_gradientn(colours = topo.colors(20)) ggsave(uhi_u,filename = paste0(wd,"UHImax/fig/uhi_u.png")) ggsave(uhi_DTR,filename = paste0(wd,"UHImax/fig/uhi_DTR.png")) ggsave(uhi_S,filename = paste0(wd,"UHImax/fig/uhi_S.png")) uhi_u uhi_DTR uhi_S caret::R2((2-uhimax_Utrecht$svf-uhimax_Utrecht$fveg)*uhimax_Utrecht$meteo,uhimax_Utrecht$UHImeasured) caret::RMSE((2-uhimax_Utrecht$svf-uhimax_Utrecht$fveg)*uhimax_Utrecht$meteo,uhimax_Utrecht$UHImeasured) R2.exp <- expression(paste(" ",R^2 ,"= 0.38")) uhimax_corr<-ggplot(uhimax_Utrecht,aes(UHImeasured,(2-svf-fveg)*meteo,colour=factor(stn))) + geom_point(size=3) + geom_abline() + xlim(0,10) + ylim(0,10) + coord_fixed() + xlab("UHImax measured") + ylab("UHImax calculated") + scale_color_manual(values=c('orange','yellow','green' )) + guides(color=guide_legend("Station ID"))+ annotate("text",x=1,y=10,label=R2.exp) + annotate("text",x=1,y=9.4,label=paste0("RMSE= ",round(wur.RMSE,2))) uhimax_corr ggsave(uhimax_corr, filename="C:/Users/marie/OneDrive/Afbeeldingen/wunderground_utrecht.png")
Using the formula of Theeuwes for Utrecht
The formula of Theeuwes(2017) was derived for:
- 0 < fveg < 0.4
- 0.2 < SVF < 0.9
Outside of this range the formula has not been tested.
svf.r<-raster::resample(grd.svf,grd.fveg,method="bilinear") values(svf.r)[values(svf.r)<0.2 | values(svf.r)>0.9] = NA values(grd.fveg)[values(grd.fveg)<0 | values(grd.fveg)>0.4] = NA city_center<-extent(c(5.1,5.13,52.08,52.098)) city_center<-raster(city_center) values(city_center)<-1 crs(city_center)<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") city_center<-projectRaster(city_center,crs=crs(grd.fveg)) st<-(2-svf.r-grd.fveg)*mean(uhimax_Utrecht$meteo) st<-crop(st,city_center) mapview(st)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.