통계 기법 지상 측정 기반 분포 지도 (농도, 오차)

draw_ground_asia_pm25.R

################(1) 패키지 및 자료 불러오기########################

install.packages('maptools');
install.packages('rgdal');
install.packages('sp');
install.packages('tmap');
install.packages('gstat');
install.packages('raster');
install.packages('rgeos');
install.packages('dplyr');
install.packages('caret');

library('maptools');
library('rgdal');
library('sp');
library('tmap');
library('gstat');
library('raster');
library('rgeos');
library('dplyr');
library('caret')

setwd("~/GroundPM2.5_Asia");

# 데이터 파일을 불러오는 코드입니다.
# ‘예시_한국/중국‘ 파일은 실제로 측정소 농도 자료가 실시간으로 들어오면 해당 파일로 넣으시면 됩니다.
# ’측정소좌표_한국/중국’ 파일은 최신 측정소 정보(위치 및 측정소명) 업데이트한이므로 제가 드린 파일로 계속 사용하시면 됩니다.
data_KR <- read.csv("예시_한국.csv", fileEncoding="CP949")
location_KR <- read.csv("측정소좌표_한국.csv", fileEncoding="CP949")
data_CN <- read.csv("예시_중국.csv", fileEncoding="CP949")
location_CN <- read.csv("측정소좌표_중국.csv", fileEncoding="CP949")

# 실시간 농도 자료+측정소 정보(좌표, 이름) 결합 코드
PM_KR <- merge(data_KR, location_KR, by="측정소명")
PM_CN <- merge(data_CN, location_CN, by="측정소명")

# 한국 자료, 중국 자료 결합 코드
PM <- rbind(PM_KR, PM_CN)
# raw data 백업용
PMraw <- PM

# 결측자료나, 가끔씩 중국 자료에서 결측값이 다르게 표시(예: NaN)으로 되어있을 때, 처리하는 코드입니다.
# NaN 부분은 상황에 맞게 수정하면 됩니다.
PM <- PM[!(PM$PM25 == "NaN"),]
PM <- na.omit(PM)
# 오차 계산에 필요한 raw data
PM_err <- PM

################(2) 지상관측기반 농도지도 작성 및 그리드별 농도 저장########################
# 분포지도 작성에 필요한 동아시아 지도 및 적용반경 shp 파일 불러오기입니다. 
# 제가 새로 보내드린 GroundPM2.5_Asia 폴더 안 GIS 폴더 안에 있는 파일 그대로 사용하시면 됩니다.
CRS_WGS <- CRS("+proj=longlat +datum=WGS84")

EastAsia <- readShapePoly("GIS/EastAsia.shp")
buffer <- readShapePoly("GIS/buf_KrCh.shp")

# 동아시아, 적용반경, 지상 관측자료의 GIS 그림을 위한 좌표설정 코드
proj4string(EastAsia) <- CRS("+proj=longlat +datum=WGS84")
proj4string(buffer) <- CRS("+proj=longlat +datum=WGS84")
coordinates(PM) <- ~lon+lat
proj4string(PM) <- CRS("+proj=longlat +datum=WGS84")

# 보간법 계산에 필요한 grid 생성 코드
grid <- as.data.frame(spsample(EastAsia, "regular", n=50000));
grid_extract <- as.data.frame(spsample(EastAsia, "regular", n=50000));
names(grid) <- c("lon", "lat");
coordinates(grid) <- c("lon", "lat");
gridded(grid) <- TRUE;
fullgrid(grid) <- TRUE;

# 좌표계 통일
proj4string(grid) <- proj4string(PM);

# IDW 보간법 적용(분포지도 작성)
data.idw <- gstat::idw(PM25~1, PM, newdata=grid, idp=2.0);
data.raster <- raster(data.idw);

# 동아시아 크기에 맞게 자른 후, 적용 반경안에 농도만 자르는 코드
PMplot <- mask(data.raster, buffer);
PMplot2 <- mask(PMplot, EastAsia);
PMconc <- extract(PMplot2, grid_extract)
Extract_PM_conc <- cbind(grid_extract, PMconc)
names(Extract_PM_conc) <- c('lat', 'lon', 'PM25conc')

# 지상관측 기반 농도분포지도의 grid 별 농도값을 csv 파일로 저장하는 코드입니다. “ ” 사이에 원하시는 파일명으로 수정가능합니다. 
write.csv(Extract_PM_conc,"Ground observ_conc_EA.csv")

################(3) 지상관측기반 오차지도 작성 및 그리드별 농도 저장########################

# 오차 계산을 위해 무작위로 지상 측정소 자료를 7:3으로 나누는 코드입니다.
div <- createDataPartition(PM_err$PM25, p=0.7, list=F)

Sample <- PM_err[div,]
Val <- PM_err[-div,]
Sample_0 <- Sample
Val_0 <- Val

# 나눠진 Sample 측정소와 Validation 측정소 자료에 대한 좌표계 설정
coordinates(Sample) <- ~lon+lat
proj4string(Sample) <- CRS("+proj=longlat +datum=WGS84")
coordinates(Val) <- ~lon+lat
proj4string(Val) <- CRS("+proj=longlat +datum=WGS84")

# 오차지도를 작성하기 위한 grid 생성 코드
grid_error <- as.data.frame(spsample(EastAsia, "regular", n=50000));
names(grid_error) <- c("lon", "lat");
coordinates(grid_error) <- c("lon", "lat");
gridded(grid_error) <- TRUE;
fullgrid(grid_error) <- TRUE;

# 좌표계 통일
proj4string(grid_error) <- proj4string(Sample);
proj4string(grid_error) <- proj4string(Val);

# Sample 측정소의 IDW 적용을 통한 분포지도 작성
Sam.idw <- gstat::idw(PM25~1, Sample, newdata=grid_error, idp=2.0);
Sam.raster <- raster(Sam.idw);

# Validation 측정소 지점에서 측정값과 Sample 측정소 자료로 예측한 예측값을 비교하여 오차(Error) 계산
Sam.pre <- extract(Sam.raster, Val)
Error <- cbind(Val_0, Sam.pre)
Ground_err <- data.frame(abs(Error$PM25-Error$Sam.pre))
names(Ground_err) <- c('PM25err')
Ground_err <- cbind(Error, Ground_err)

# 오차 자료 백업
Ground_err0 <- Ground_err
# 계산된 오차 농도 자료에 좌표계 설정
coordinates(Ground_err) <- ~lon+lat
proj4string(Ground_err) <- CRS("+proj=longlat +datum=WGS84")
proj4string(grid_error) <- proj4string(Ground_err);
# 오차 농도로 오차지도 작성
error.idw <- gstat::idw(PM25err~1, Ground_err, newdata=grid_error, idp=2.0);
error.raster <- raster(error.idw);

# 오차지도에 동아시아 크기로 자르고, 적용 반경안에 농도만 자르는 코드
PMerr <- mask(error.raster, buffer);
PMerr2 <- mask(PMerr, EastAsia);
PMerror <- extract(PMerr2, grid_extract)
Extract_PM_err <- cbind(grid_extract, PMerror)
names(Extract_PM_err) <- c('lat', 'lon', 'PM25err')

# 지상관측 오차 분포지도의 grid 별 농도값을 csv 파일로 저장하는 코드입니다. “ ” 사이에 원하시는 파일명으로 수정가능합니다. 
write.csv(Extract_PM_err,"Ground observ_error_EA.csv")

################(4) 지상관측기반 농도 및 오차지도 저장########################

# 농도분포지도 그리기
tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_shape(PMplot2)+tm_raster(n=30, style="fixed", breaks=c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 90, 100, 110, 120, 500), 
                                                                    palette=c("mediumblue", "blue", "deepskyblue2", "cyan", "green", "greenyellow", "yellow", "gold", "orange2", "darkorange2",  "red", "red2", "red3", "orangered4", "red4", "coral4", "black"), legend.reverse=TRUE, 
                                                                    title = "PM2.5 conc. (ug/m3)")+tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_legend(legend.outside=TRUE, legend.title.size=2, legend.text.size=1.5);
# 오차분포지도 그리기
tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_shape(PMerr2)+tm_raster(n=30, style="fixed", breaks=c(0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 500), 
                                                                   palette=c("mediumblue", "blue", "deepskyblue2", "cyan", "green", "greenyellow", "yellow", "gold", "orange2", "darkorange2",  "red", "red2", "red3", "orangered4", "red4", "coral4", "black"), legend.reverse=TRUE, 
                                                                   title = "PM2.5 error. (ug/m3)")+tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_legend(legend.outside=TRUE, legend.title.size=2, legend.text.size=1.5);

PMmap_conc <- tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_shape(PMplot2)+tm_raster(n=30, style="fixed", breaks=c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 90, 100, 110, 120, 500), 
                                                                                  palette=c("mediumblue", "blue", "deepskyblue2", "cyan", "green", "greenyellow", "yellow", "gold", "orange2", "darkorange2",  "red", "red2", "red3", "orangered4", "red4", "coral4", "black"), legend.reverse=TRUE, 
                                                                                  title = "PM2.5 conc. (ug/m3)")+tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_legend(legend.outside=TRUE, legend.title.size=2, legend.text.size=1.5);

PMmap_error <- tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_shape(PMerr2)+tm_raster(n=30, style="fixed", breaks=c(0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 500), 
                                                                                  palette=c("mediumblue", "blue", "deepskyblue2", "cyan", "green", "greenyellow", "yellow", "gold", "orange2", "darkorange2",  "red", "red2", "red3", "orangered4", "red4", "coral4", "black"), legend.reverse=TRUE, 
                                                                                  title = "PM2.5 error. (ug/m3)")+tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_legend(legend.outside=TRUE, legend.title.size=2, legend.text.size=1.5);

tmap_save(PMmap_conc, width=20, height=20, units="cm", dpi=300, filename="PMmap_conc.png")	
tmap_save(PMmap_error, width=20, height=20, units="cm", dpi=300, filename="PMmap_error.png")	

################(1) 패키지 및 자료 불러오기########################

install.packages('maptools');
install.packages('rgdal');
install.packages('sp');
install.packages('tmap');
install.packages('gstat');
install.packages('raster');
install.packages('rgeos');
install.packages('dplyr');
install.packages('caret');
install.packages('sf');

library('maptools');
library('rgdal');
library('sp');
library('tmap');
library('gstat');
library('raster');
library('rgeos');
library('dplyr');
library('caret')

setwd("/app/ground-asia-pm25-conc-err");
Sys.umask("000")

# 데이터 파일 sample.
# data_KR <- read.csv("data-korea.csv")
# data_CN <- read.csv("data-china.csv")

# 어제와 오늘 두개의 파일을 처리한다.
china_files <- list.files('/finepart/ground-china', recursive = TRUE)
china_last_day_file <- china_files[length(china_files)]
china_last_1_day_file <- china_files[length(china_files)-1]
china_last_file <- paste("/finepart/ground-china", china_last_day_file, sep='/')
china_last_1_file <- paste("/finepart/ground-china", china_last_1_day_file, sep='/')

china_day_data <- read.csv(china_last_file)
china_day_1_data <- read.csv(china_last_1_file)
china_day_data2 <- subset(china_day_data, type == "PM2.5", select = -c(type))
china_day_1_data2 <- subset(china_day_1_data, type == "PM2.5", select = -c(type))
# 두 파일을 수직으로 붙인다.
china_data2 <- rbind(china_day_1_data2, china_day_data2)
# 아래 두 줄만 남긴다.
china_data3 <- tail(china_data2, n = 2)
# 첫째 줄이 1시간 전, 둘째 줄이 마지막 시간.
china_first <- head(china_data3, n = 1)
china_second <- tail(china_data3, n = 1)

# 1시간 전은 무시하고 마지막 시간만 해보자.
china4 <- t(china_second)
china4_date_part <- head(china4, n = 2)

# df --> string
china_date <- paste(as.character(china4_date_part), collapse = '')
# remove date part
china5 <- subset(china_second, select = -c(date, hour))
china6 <- data.frame(t(china5))
# index to first column
data_CN <- cbind(loc = rownames(china6), data.frame(china6, row.names=NULL))
# rename column names
colnames(data_CN) <- c("측정소명", "PM25")

# 중국 데이터 준비 끝.

# 이제 한국. 시차가 1시간이므로 중국이 01시 이면 우리는 02시 파일을 선택해야 함.
china_time <- strptime(china_date, format = '%Y%m%d%H')
korea_time <- china_time + 60*60
korea_date <- format(korea_time, '%Y%m%d%H')
korea_file <- paste('/finepart/ground-korea', substring(korea_date,0,4), substring(korea_date,0,6), paste(korea_date, "-latlon.csv", sep=''), sep='/')

file.exists(korea_file)
if (file.exists(korea_file)) {
  data_KR_1 <- read.csv(korea_file)
  data_KR <- data_KR_1[, c('지역', 'pm25Value')]
  colnames(data_KR) <- c("측정소명", "PM25")
  
  # 측정소 위치 자료. 2022.7월에 업데이트된 위치.
  location_KR <- read.csv("location-korea.csv")
  location_CN_1 <- read.csv("location-china.csv")
  
  # 데이터의 측정소명은 앞에 'X'가 붙어 있다. 그래서 위치 정보의 측정소명에도 'X'를 붙여 준다.
  location_CN_1$Xloc <-paste('X', location_CN_1$측정소명, sep = '')
  location_CN <- location_CN_1[, c('Xloc', 'lon', 'lat')]
  colnames(location_CN) <- c("측정소명", 'lon', 'lat')
  
  # 실시간 농도 자료+측정소 정보(좌표, 이름) 결합 코드
  PM_KR <- merge(data_KR, location_KR, by="측정소명")
  PM_CN <- merge(data_CN, location_CN, by="측정소명")
  
  # 한국 자료, 중국 자료 결합 코드
  PM <- rbind(PM_KR, PM_CN)
  # PM25 칼럼을 numeric으로 변경
  PM$PM25 <- as.numeric(as.character(PM$PM25))
  # raw data 백업용
  PMraw <- PM
  
  # 결측자료나, 가끔씩 중국 자료에서 결측값이 다르게 표시(예: NaN)으로 되어있을 때, 처리하는 코드입니다.
  # NaN 부분은 상황에 맞게 수정하면 됩니다.
  PM <- PM[!(PM$PM25 == "NaN"),]
  PM <- na.omit(PM)
  # 오차 계산에 필요한 raw data
  PM_err <- PM
  
  ################(2) 지상관측기반 농도지도 작성 및 그리드별 농도 저장########################
  # 분포지도 작성에 필요한 동아시아 지도 및 적용반경 shp 파일 불러오기입니다. 
  # 제가 새로 보내드린 GroundPM2.5_Asia 폴더 안 GIS 폴더 안에 있는 파일 그대로 사용하시면 됩니다.
  CRS_WGS <- CRS("+proj=longlat +datum=WGS84")
  
  EastAsia <- readShapePoly("GIS/EastAsia.shp")
  buffer <- readShapePoly("GIS/Buf_KrCh.shp")
  
  # 동아시아, 적용반경, 지상 관측자료의 GIS 그림을 위한 좌표설정 코드
  proj4string(EastAsia) <- CRS("+proj=longlat +datum=WGS84")
  proj4string(buffer) <- CRS("+proj=longlat +datum=WGS84")
  coordinates(PM) <- ~lon+lat
  proj4string(PM) <- CRS("+proj=longlat +datum=WGS84")
  
  # 보간법 계산에 필요한 grid 생성 코드
  grid <- as.data.frame(spsample(EastAsia, "regular", n=50000));
  grid_extract <- as.data.frame(spsample(EastAsia, "regular", n=50000));
  names(grid) <- c("lon", "lat");
  coordinates(grid) <- c("lon", "lat");
  gridded(grid) <- TRUE;
  fullgrid(grid) <- TRUE;
  
  # 좌표계 통일
  proj4string(grid) <- proj4string(PM);
  
  # IDW 보간법 적용(분포지도 작성)
  data.idw <- gstat::idw(PM25~1, PM, newdata=grid, idp=2.0);
  data.raster <- raster(data.idw);
  
  # 동아시아 크기에 맞게 자른 후, 적용 반경안에 농도만 자르는 코드
  PMplot <- mask(data.raster, buffer);
  PMplot2 <- mask(PMplot, EastAsia);
  PMconc <- extract(PMplot2, grid_extract)
  Extract_PM_conc <- cbind(grid_extract, PMconc)
  names(Extract_PM_conc) <- c('lat', 'lon', 'PM25conc')

  dir.create(paste('/finepart/ground-asia-pm25-conc-err/', substr(korea_date, 0, 4), substr(korea_date, 0, 6), sep = '/'), showWarnings = FALSE, recursive = TRUE);  
  # 지상관측 기반 농도분포지도의 grid 별 농도값을 csv 파일로 저장하는 코드입니다. “ ” 사이에 원하시는 파일명으로 수정가능합니다. 
  conc_csv_file = paste('/finepart/ground-asia-pm25-conc-err/', 
                        substr(korea_date, 0, 4), '/', substr(korea_date, 0, 6), 
                        '/ground-asia-pm25-conc-', korea_date, ".csv", sep = '')
  write.csv(Extract_PM_conc, conc_csv_file)
  
  ################(3) 지상관측기반 오차지도 작성 및 그리드별 농도 저장########################
  
  # 오차 계산을 위해 무작위로 지상 측정소 자료를 7:3으로 나누는 코드입니다.
  div <- createDataPartition(PM_err$PM25, p=0.7, list=F)
  
  Sample <- PM_err[div,]
  Val <- PM_err[-div,]
  Sample_0 <- Sample
  Val_0 <- Val
  
  # 나눠진 Sample 측정소와 Validation 측정소 자료에 대한 좌표계 설정
  coordinates(Sample) <- ~lon+lat
  proj4string(Sample) <- CRS("+proj=longlat +datum=WGS84")
  coordinates(Val) <- ~lon+lat
  proj4string(Val) <- CRS("+proj=longlat +datum=WGS84")
  
  # 오차지도를 작성하기 위한 grid 생성 코드
  grid_error <- as.data.frame(spsample(EastAsia, "regular", n=50000));
  names(grid_error) <- c("lon", "lat");
  coordinates(grid_error) <- c("lon", "lat");
  gridded(grid_error) <- TRUE;
  fullgrid(grid_error) <- TRUE;
  
  # 좌표계 통일
  proj4string(grid_error) <- proj4string(Sample);
  proj4string(grid_error) <- proj4string(Val);
  
  # Sample 측정소의 IDW 적용을 통한 분포지도 작성
  Sam.idw <- gstat::idw(PM25~1, Sample, newdata=grid_error, idp=2.0);
  Sam.raster <- raster(Sam.idw);
  
  # Validation 측정소 지점에서 측정값과 Sample 측정소 자료로 예측한 예측값을 비교하여 오차(Error) 계산
  Sam.pre <- extract(Sam.raster, Val)
  Error <- cbind(Val_0, Sam.pre)
  Ground_err <- data.frame(abs(Error$PM25-Error$Sam.pre))
  names(Ground_err) <- c('PM25err')
  Ground_err <- cbind(Error, Ground_err)
  
  # 오차 자료 백업
  Ground_err0 <- Ground_err
  # 계산된 오차 농도 자료에 좌표계 설정
  coordinates(Ground_err) <- ~lon+lat
  proj4string(Ground_err) <- CRS("+proj=longlat +datum=WGS84")
  proj4string(grid_error) <- proj4string(Ground_err);
  # 오차 농도로 오차지도 작성
  error.idw <- gstat::idw(PM25err~1, Ground_err, newdata=grid_error, idp=2.0);
  error.raster <- raster(error.idw);
  
  # 오차지도에 동아시아 크기로 자르고, 적용 반경안에 농도만 자르는 코드
  PMerr <- mask(error.raster, buffer);
  PMerr2 <- mask(PMerr, EastAsia);
  PMerror <- extract(PMerr2, grid_extract)
  Extract_PM_err <- cbind(grid_extract, PMerror)
  names(Extract_PM_err) <- c('lat', 'lon', 'PM25err')
  
  # 지상관측 오차 분포지도의 grid 별 농도값을 csv 파일로 저장하는 코드입니다. “ ” 사이에 원하시는 파일명으로 수정가능합니다. 
  err_csv_file = paste('/finepart/ground-asia-pm25-conc-err/', 
                        substr(korea_date, 0, 4), '/', substr(korea_date, 0, 6), 
                        '/ground-asia-pm25-err-', korea_date, ".csv", sep = '')
  write.csv(Extract_PM_err, err_csv_file)
  
  ################(4) 지상관측기반 농도 및 오차지도 저장########################
  
  # 농도분포지도 그리기
  tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_shape(PMplot2)+tm_raster(n=30, style="fixed", breaks=c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 90, 100, 110, 120, 500), 
                                                                      palette=c("mediumblue", "blue", "deepskyblue2", "cyan", "green", "greenyellow", "yellow", "gold", "orange2", "darkorange2",  "red", "red2", "red3", "orangered4", "red4", "coral4", "black"), legend.reverse=TRUE, 
                                                                      title = "PM2.5 conc. (ug/m3)")+tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_legend(legend.outside=TRUE, legend.title.size=2, legend.text.size=1.5);
  # 오차분포지도 그리기
  tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_shape(PMerr2)+tm_raster(n=30, style="fixed", breaks=c(0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 500), 
                                                                     palette=c("mediumblue", "blue", "deepskyblue2", "cyan", "green", "greenyellow", "yellow", "gold", "orange2", "darkorange2",  "red", "red2", "red3", "orangered4", "red4", "coral4", "black"), legend.reverse=TRUE, 
                                                                     title = "PM2.5 error. (ug/m3)")+tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_legend(legend.outside=TRUE, legend.title.size=2, legend.text.size=1.5);
  
  PMmap_conc <- tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_shape(PMplot2)+tm_raster(n=30, style="fixed", breaks=c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 90, 100, 110, 120, 500), 
                                                                                    palette=c("mediumblue", "blue", "deepskyblue2", "cyan", "green", "greenyellow", "yellow", "gold", "orange2", "darkorange2",  "red", "red2", "red3", "orangered4", "red4", "coral4", "black"), legend.reverse=TRUE, 
                                                                                    title = "PM2.5 conc. (ug/m3)")+tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_legend(legend.outside=TRUE, legend.title.size=2, legend.text.size=1.5);
  
  PMmap_error <- tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_shape(PMerr2)+tm_raster(n=30, style="fixed", breaks=c(0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 500), 
                                                                                    palette=c("mediumblue", "blue", "deepskyblue2", "cyan", "green", "greenyellow", "yellow", "gold", "orange2", "darkorange2",  "red", "red2", "red3", "orangered4", "red4", "coral4", "black"), legend.reverse=TRUE, 
                                                                                    title = "PM2.5 error. (ug/m3)")+tm_shape(EastAsia)+tm_polygons(alpha=0)+tm_legend(legend.outside=TRUE, legend.title.size=2, legend.text.size=1.5);
  
  tmap_save(PMmap_conc, width=20, height=20, units="cm", dpi=300, filename=paste('/finepart/ground-asia-pm25-conc-err/', substr(korea_date, 0, 4), '/', substr(korea_date, 0, 6), '/ground-asia-pm25-conc-', korea_date, ".png", sep = ''));
  tmap_save(PMmap_error, width=20, height=20, units="cm", dpi=300, filename=paste('/finepart/ground-asia-pm25-conc-err/', substr(korea_date, 0, 4), '/', substr(korea_date, 0, 6), '/ground-asia-pm25-err-', korea_date, ".png", sep = ''));
} else {
  # no korea file
}
korea_date