GCMC_RScripts/Create_AP_PMComponent_measures.R at main · WillhKessler/GCMC_RScripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
library('terra')
library('tools')


projection_info<-"EPSG:4326"

files<-list.files("S:\\GCMC\\tmp\\pmcomp_temp",pattern=".*.rds$",recursive=TRUE,full.names=TRUE)
done<-list.files("S:\\GCMC\\tmp\\pmcomp_temp",pattern=".*.tif$",recursive=TRUE,full.names=TRUE)
files<-files[!basename(file_path_sans_ext(files)) %in% basename(file_path_sans_ext(done))]
grids<- list.files("S:\\GCMC\\tmp\\pmcomp_temp\\grids",pattern="*.gpkg$",recursive=TRUE,full.names=TRUE)

for(g in grids){
  gridname<- file_path_sans_ext(basename(g))
  region_files<-files[grep(gridname,files)]
  for(compfile in region_files){
    dat<-readRDS(compfile)
    #dat<-as.data.frame(dat)

    datvect<-vect(g)
    datvect$value<-dat

    datvect<-terra::project(datvect,"ESRI:102010")
    if(grepl("non-urban",gridname)){
      datrast<-rast(extend(ext(datvect),500),res=c(1000,1000),crs="ESRI:102010")
      }else{
    datrast<-rast(extend(ext(datvect),25),res=c(50,50),crs="ESRI:102010")
    }

    datrast<-rasterize(x=datvect,
                       y=datrast,
                       field="value",
                       fun=mean,
                       filename= paste0(dirname(dirname(file_path_sans_ext(compfile))),"\\",basename(file_path_sans_ext(compfile)),".tif"),overwrite=TRUE)
  }
}

######Confirm all rasters are valid:

files<-list.files("S:\\GCMC\\tmp\\pmcomp_temp",pattern=".*.rds$",recursive=TRUE,full.names=TRUE)
files[basename(files) %in% files[duplicated(basename(files))]]


######Create single layer coverages

library(terra)
require(tools)

files<-list.files("S:\\GCMC\\tmp\\pmcomp_temp",recursive=TRUE,pattern="*.tif$",full.names = TRUE)

patterns = unlist(unique(lapply(X = strsplit(file_path_sans_ext(basename(files)),split="_"),
                                FUN = function(x){paste0(unlist(x)[1:4],collapse = "_")})))


for(i in patterns){
  urbanfiles<-list.files("S:/GCMC/Data/AirPollution/PM25_Components/urban",
                         recursive=TRUE,
                         pattern=paste0(i,".*.tif$"),
                         full.names = TRUE)
  nonurbanfiles<-list.files("S:/GCMC/Data/AirPollution/PM25_Components/non_urban",
                            recursive=TRUE,
                            pattern=paste0(i,".*.tif$"),
                            full.names = TRUE)


  # Read in files as rasters, one for urban, one for non-urban
  urbanrast<-rast(urbanfiles[1])
  nonurbanrast<-rast(nonurbanfiles[1])

  # extend both rasters so their extents match
  #urbanrast<- extend(urbanrast,nonurbanrast,snap = "near")
  #nonurbanrast<-extend(nonurbanrast,urbanrast,snap = "near")

  # create an empty raster of the non urban extent and resolution
  nonurbanrastnofill<-nonurbanrast
  values(nonurbanrastnofill)<-NA

  # Resample the urban rast to the non-urban resolution
  # resample_urbanrast<-resample(urbanrast,nonurbanrast,method="cubicspline",threads=TRUE,filename="S:/GCMC/tmp/resample.tif")
  resample_urbanrast<-resample(urbanrast,nonurbanrast,method="cubicspline",threads=TRUE)

  # Mosaic the urban and non-urban rasters together
  # PM25_comp<-mosaic(nonurbanrast,resample_urbanrast,fun="first",filename="S:/GCMC/tmp/mosaic.tif")
  PM25_comp<-mosaic(nonurbanrast,resample_urbanrast,fun="first")

  # Use a void fill to fill in any remaining NA values using an average of surrounding values
  PM25_comp<-focal(PM25_comp,w=3,fun=mean,na.policy="only",na.rm=T,filename=file.path(dirname(dirname(dirname(dirname(urbanfiles[1])))),paste0(paste0(unlist(strsplit(i[[1]],"_"))[-1],collapse="_"),".tif")),
                   overwrite=TRUE)

  #outrast<- rast(list(urbanrast,nonurbanrast,PM25_comp))
  #writeRaster(outrast,filename=file.path(dirname(dirname(dirname(dirname(urbanfiles[1])))),paste0(paste0(unlist(strsplit(i[[1]],"_"))[-1],collapse="_"),".tif")),
  #            overwrite=TRUE))


}