finep hybrid 생산

draw_hybrid.R

# 1. 세팅 단계
# 1-1. 필요한 패키지 불러오기
library('maptools');
library('rgdal');
library('sp');
library('tmap');
library('gstat');
library('raster');
library('rgeos');

# 1-2. 경로 설정
setwd("/app/PM2.5map_Korea");
Sys.umask("000")

# 1-3. 대한민국 지도 불러온 후 좌표계 설정
CRS_WGS <- CRS("+proj=longlat +datum=WGS84")
KoreaTM <- readShapePoly("GIS/temp.shp");	 ############ Input1: 한반도 GIS shp 파일 ###################

proj4string(KoreaTM) <- CRS("+proj=tmerc +lat_0=38 +lon_0=127 +k=1 +x_0=200000 +y_0=500000 +ellps=bessel +units=m")
KoreaWGS <- spTransform(KoreaTM, CRS_WGS)

# 1-4. 위성 측정 좌표(동아시아) 좌표계 설정
Sat_grid <- read.csv("latlon_goci.csv")
Sat_grid0 <- Sat_grid
coordinates(Sat_grid)=~lon+lat
Sat_grid <- SpatialPoints(Sat_grid, proj4string=CRS_WGS)

# 2. 위성 관측지도 작성(대한민국)
# 2-1. 위성 측정값 불러온 후 좌표계 설정(대한민국 범위만 자르기)
Sat_conc_raw <- read.csv("/finepart/ground-korea/{{year-dir}}/{{month-dir}}/{{conc-csv-filename}}");	 ############ Input2: 위성예측농도 csv파일 ###################
Sat_conc <- cbind(Sat_grid0, Sat_conc_raw)
names(Sat_conc) <- c( 'lat', 'lon', 'Sat_PM2.5_conc')
Sat_conc0 <- Sat_conc
coordinates(Sat_conc)=~lon+lat
Sat_conc <- SpatialPoints(Sat_conc, proj4string=CRS_WGS)
Sat_conc_Korea <- cbind(Sat_conc0, extract(KoreaWGS, Sat_conc))
Sat_conc_Korea <- na.omit(Sat_conc_Korea)

location <- Sat_conc_Korea
location0 <- location[1:2]
coordinates(location)=~lon+lat
location <- SpatialPoints(location, proj4string=CRS_WGS)

# 2-2. 위성관측값 지도 작성
# 2-2-1. 결측값 처리
if(length(which(Sat_conc_Korea$Sat_PM2.5_conc==-9999))==length(Sat_conc_Korea$lat)) {
    Sat_conc_Korea <- replace(Sat_conc_Korea$Sat_PM2.5_conc, Sat_conc_Korea$Sat_PM2.5_conc==-9999,0)
}

Sat_conc_Korea <- replace(Sat_conc_Korea, Sat_conc_Korea==-9999,NA)
Sat_conc_Korea <- Sat_conc_Korea[1:3]
Sat_conc_Korea <- na.omit(Sat_conc_Korea)

#2-2-2. 결측값 처리 된 위성 지도 작성 및 농도 추출
coordinates(Sat_conc_Korea)=~lon+lat
proj4string(Sat_conc_Korea) <- CRS("+proj=longlat +datum=WGS84")

grid1 <- as.data.frame(spsample(location, "regular", n=50000));
names(grid1) <- c("lon", "lat");
coordinates(grid1) <- c("lon", "lat");
gridded(grid1) <- TRUE;
fullgrid(grid1) <- TRUE;

proj4string(grid1) <- proj4string(Sat_conc_Korea);

Sat_conc_Korea.idw <- gstat::idw(Sat_PM2.5_conc~1, Sat_conc_Korea, newdata=grid1, idp=2.0);
Sat_conc_Korea.raster <- raster(Sat_conc_Korea.idw);

SatPMmap <- mask(Sat_conc_Korea.raster, KoreaWGS)
Sat_conc.level <- extract(Sat_conc_Korea.raster, location)

Sat_concPM2.5 <- cbind(location0, Sat_conc.level)

# 3. 위성 오차값 지도
# 3-1. 위성 오차값 불러온 후 좌표계 설정(대한민국 범위만 자르기)
Sat_err_raw <- read.csv("/finepart/ground-korea/{{year-dir}}/{{month-dir}}/{{err-csv-filename}}")	############ Input3: 위성오차값 csv파일###################
Sat_err <- cbind(Sat_err_raw[8:9], Sat_err_raw[2])
names(Sat_err) <- c('lat', 'lon', 'Sat_PM2.5_err')
Sat_err0 <- Sat_err
coordinates(Sat_err)=~lon+lat
Sat_err <- SpatialPoints(Sat_err, proj4string=CRS_WGS)
Sat_err_Korea <- cbind(Sat_err0, extract(KoreaWGS, Sat_err))
Sat_err_Korea <- na.omit(Sat_err_Korea)

# 3-2. 위성오차값 지도 작성
# 3-2-1. 결측 값 처리
if(length(which(Sat_err_Korea$Sat_PM2.5_err==-9999))==length(Sat_err_Korea$lat)) {
    Sat_err_Korea <- replace(Sat_err_Korea$Sat_PM2.5_err, Sat_err_Korea$Sat_PM2.5_err==-9999, 0)
}

Sat_err_Korea <- replace(Sat_err_Korea, Sat_err_Korea==-9999,NA)
Sat_err_Korea <- na.omit(Sat_err_Korea)
Sat_err_Korea <- cbind(Sat_err_Korea[1:2], abs(Sat_err_Korea[3]))

#3-2-2. 결측값 처리 된 위성 오차 지도 작성
coordinates(Sat_err_Korea)=~lon+lat
proj4string(Sat_err_Korea) <- CRS("+proj=longlat +datum=WGS84")

grid2 <- as.data.frame(spsample(location, "regular", n=50000));
names(grid2) <- c("lon", "lat");
coordinates(grid2) <- c("lon", "lat");
gridded(grid2) <- TRUE;
fullgrid(grid2) <- TRUE;

proj4string(grid2) <- proj4string(Sat_err_Korea);

Sat_err_Korea.idw <- gstat::idw(Sat_PM2.5_err~1, Sat_err_Korea, newdata=grid2, idp=2.0);
Sat_err_Korea.raster <- raster(Sat_err_Korea.idw);

SatPMerrmap <- mask(Sat_err_Korea.raster, KoreaWGS)

#3-2-3. 위성 오차값 추출(격자별)
Sat_err.level <- extract(Sat_err_Korea.raster, location)

Sat_errPM2.5 <- cbind(location0, Sat_err.level)



# 4. 지상 관측지도 작성(대한민국)
# 4-1. 대한민국 지상측정소 자료 불러온 후 좌표계 설정
Ground_conc <- cbind(Sat_err_raw[8:9], Sat_err_raw[6])
names(Ground_conc) <- c('lat', 'lon', 'Ground_PM2.5_conc')
Ground_conc0 <- Ground_conc
coordinates(Ground_conc)=~lon+lat
Ground_conc <- SpatialPoints(Ground_conc, proj4string=CRS_WGS)
Ground_conc_Korea <- cbind(Ground_conc0, extract(KoreaWGS, Ground_conc))
Ground_conc_Korea <- na.omit(Ground_conc_Korea)

# 4-2. 지상관측값 지도 작성
# 4-2-1. 결측값 처리
if(length(which(Ground_conc_Korea$Ground_PM2.5_err==-9999))==length(Ground_conc_Korea$lat)) {
    Ground_conc_Korea <- replace(Ground_conc_Korea$Ground_PM2.5, Ground_conc_Korea$Ground_PM2.5==-9999,0)
}

Ground_conc_Korea <- replace(Ground_conc_Korea, Ground_conc_Korea==-9999, NA)
Ground_conc_Korea <- na.omit(Ground_conc_Korea)
Ground_conc_Korea0 <- Ground_conc_Korea[1:3]

#4-2-2. 결측값 처리 된 지상 관측값 지도 작성
coordinates(Ground_conc_Korea)=~lon+lat
proj4string(Ground_conc_Korea) <- CRS("+proj=longlat +datum=WGS84")

grid3 <- as.data.frame(spsample(location, "regular", n=50000));
names(grid3) <- c("lon", "lat");
coordinates(grid3) <- c("lon", "lat");
gridded(grid3) <- TRUE;
fullgrid(grid3) <- TRUE;

proj4string(grid3) <- proj4string(Ground_conc_Korea);

Ground_conc_Korea.idw <- gstat::idw(Ground_PM2.5_conc~1, Ground_conc_Korea, newdata=grid3, idp=2.0);
Ground_conc_Korea.raster <- raster(Ground_conc_Korea.idw);

GroundPMmap <- mask(Ground_conc_Korea.raster, KoreaWGS)

#4-2-3. 지상 관측값 저장(지점별)
Ground.conc.level <- extract(Ground_conc_Korea.raster, location)

Ground_concPM2.5 <- cbind(location0, Ground.conc.level)
write.csv(Ground_concPM2.5, "/finepart/hybrid/{{year-dir}}/{{month-dir}}/Ground_conc_Korea_{{datetime}}.csv")	############ Output1: 지상예측농도(격자별) csv 파일 ###################

# 5. 지상 오차지도 작성(대한민국)
# 5-1. 지점별 예측값 계산
Ground_length <- length(Ground_conc_Korea0$lat)

var <- c(1:3)
pre <- as.numeric(var)

for(i in 1:Ground_length)
{	point <- Ground_conc_Korea0[i,]
	point0 <- point[1:2]
	Ground_pre <- Ground_conc_Korea0[-c(i),]

	coordinates(Ground_pre)=~lon+lat;
	coordinates(point)=~lon+lat;

	proj4string(Ground_pre) <- CRS("+proj=longlat +datum=WGS84")
	proj4string(point) <- CRS("+proj=longlat +datum=WGS84")

	grid4 <- as.data.frame(spsample(location, "regular", n=50000));
	names(grid4) <- c("lon", "lat");
	coordinates(grid4) <- c("lon", "lat");
	gridded(grid4) <- TRUE;
	fullgrid(grid4) <- TRUE;

	proj4string(grid4) <- proj4string(Ground_pre);

	Ground_pre.idw <- gstat::idw(Ground_PM2.5_conc~1, Ground_pre, newdata=grid4, idp=2.0);
	Ground_pre.raster <- raster(Ground_pre.idw);

	pre.level <- extract(Ground_pre.raster, point)
	pre.level2 <- cbind(point0, pre.level)
	pre <- rbind(pre, pre.level2)
}

# 5-2. 지점별 오차값 계산(실측값-예측값)
pre2 <- data.frame(pre)
pre3 <- pre2[-1,]
Ground_pre1 <- data.frame(pre3)
Ground_pre2 <- cbind(Ground_pre1, Ground_conc_Korea0[3])
Ground_pre2 <- na.omit(Ground_pre2)
Ground_err0 <- abs(Ground_pre2[4]-Ground_pre2[3])
names(Ground_err0) <- c('Ground_err')
Ground_err <- cbind(Ground_pre2[1:2], Ground_err0)

# 5-3. 지상 오차 지도 작성
Ground_err_map <- Ground_err
coordinates(Ground_err_map)=~lon+lat
proj4string(Ground_err_map) <- CRS("+proj=longlat +datum=WGS84")

grid5 <- as.data.frame(spsample(location, "regular", n=50000));
names(grid5) <- c("lon", "lat");
coordinates(grid5) <- c("lon", "lat");
gridded(grid5) <- TRUE;
fullgrid(grid5) <- TRUE;

proj4string(grid5) <- proj4string(Ground_err_map);

Ground_err_map.idw <- gstat::idw(Ground_err~1, Ground_err_map, newdata=grid5, idp=2.0);
Ground_err_map.raster <- raster(Ground_err_map.idw);

Grounderrmap <- mask(Ground_err_map.raster, KoreaWGS)

# 5-4. 지상 오차값 저장(지점별)
Ground_err.level <- extract(Ground_err_map.raster, location)
Ground_errPM2.5 <- cbind(location0, Ground_err.level)
write.csv(Ground_errPM2.5, "/finepart/hybrid/{{year-dir}}/{{month-dir}}/Ground_err_Korea_{{datetime}}.csv")	############ Output2: 지상오차농도(격자별) csv 파일 ###################

# 6. 통합 오차 지도 작성
# 6-1. 통합 오차 계산
Etot <- (Ground_errPM2.5[3]+Sat_errPM2.5[3])/2
Etot <- cbind(location0, Etot)
names(Etot) <- c('lat', 'lon', 'Totalerr')
Etot <- na.omit(Etot)

coordinates(Etot)=~lon+lat
proj4string(Etot) <- CRS("+proj=longlat +datum=WGS84")

grid6 <- as.data.frame(spsample(location, "regular", n=50000));
names(grid6) <- c("lon", "lat");
coordinates(grid6) <- c("lon", "lat");
gridded(grid6) <- TRUE;
fullgrid(grid6) <- TRUE;

proj4string(grid6) <- proj4string(Etot);

Etot.idw <- gstat::idw(Totalerr~1, Etot, newdata=grid6, idp=2.0);
Etot.raster <- raster(Etot.idw);

Etot.level <- extract(Etot.raster, location)
Total_err <- cbind(location0, Etot.level)
write.csv(Total_err, "/finepart/hybrid/{{year-dir}}/{{month-dir}}/Total_err_{{datetime}}.csv")	############ Output3: 통합오차농도(격자별) csv 파일 ###################

TotalerrMap <- mask(Etot.raster, KoreaWGS)
TotalerrMap2 <- tm_shape(KoreaWGS)+tm_polygons(alpha=0)+tm_shape(TotalerrMap)+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, 150),
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 err. (total) (ug/m3)")+tm_shape(KoreaWGS)+tm_polygons(alpha=0)+tm_legend(legend.outside=FALSE);
tmap_save(TotalerrMap2, width=20, height=20, units="cm", dpi=300, filename="/finepart/hybrid/{{year-dir}}/{{month-dir}}/PMtotal_err_{{datetime}}.png") 	############ Output4: 통합오차농도지도 png파일 ###################


# 7. 위성 지상 연계 지도 작성
# 7-1. 위성 지상 통합 농도 계산
PM2.5idw <- Ground_concPM2.5[3]
PM2.5sat <- Sat_concPM2.5[3]
Eidw <- Ground_errPM2.5[3]
Esat <- Sat_errPM2.5[3]

	if(length(which(Sat_concPM2.5$Sat_PM2.5_conc==0))==length(Sat_concPM2.5$lat))
	{	PM2.5idw <- GroundPM2.5[3]
		PM2.5sat <- GroundPM2.5[3]
		Eidw <- Ground_errPM2.5[3]
		Esat <- Ground_errPM2.5[3]
	}

PM2.5combined <- ((PM2.5idw*(Eidw)^-1)+(PM2.5sat*(Esat)^-1))/(((Eidw)^-1)+((Esat)^-1))

PM2.5combined2 <- cbind(location0, PM2.5combined)
names(PM2.5combined2) <- c('lat', 'lon', 'PM2.5combined')
PM2.5combined2 <- na.omit(PM2.5combined2)

coordinates(PM2.5combined2)=~lon+lat
proj4string(PM2.5combined2) <- CRS("+proj=longlat +datum=WGS84")

grid7 <- as.data.frame(spsample(location, "regular", n=50000));
names(grid7) <- c("lon", "lat");
coordinates(grid7) <- c("lon", "lat");
gridded(grid7) <- TRUE;
fullgrid(grid7) <- TRUE;

proj4string(grid7) <- proj4string(PM2.5combined2);

PM2.5combined.idw <- gstat::idw(PM2.5combined~1, PM2.5combined2, newdata=grid7, idp=2.0);
PM2.5combined.raster <- raster(PM2.5combined.idw);

PM2.5combined.level <- extract(PM2.5combined.raster, location)

PM2.5combined_data <- cbind(location0, PM2.5combined.level)
write.csv(PM2.5combined_data, "/finepart/hybrid/{{year-dir}}/{{month-dir}}/Hybrid_conc_Korea_{{datetime}}.csv")		############ Output5: 지상위성연계hybrid(격자별) csv 파일 ###################


PM2.5combinedmap <- mask(PM2.5combined.raster, KoreaWGS)
PM2.5combinedmap2 <- tm_shape(KoreaWGS)+tm_polygons(alpha=0)+tm_shape(PM2.5combinedmap)+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, 150),
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. (hybrid) (ug/m3)")+tm_shape(KoreaWGS)+tm_polygons(alpha=0)+tm_legend(legend.outside=FALSE);
tmap_save(PM2.5combinedmap2, width=20, height=20, units="cm", dpi=300, filename="/finepart/hybrid/{{year-dir}}/{{month-dir}}/PMconc_hybrid_{{datetime}}.png")		############ Output6: 지상위성연계hybrid 농도지도 png파일