매트랩을 이용해서 구글어스 위에다 그림을 그려보자.

 목표: 서울에 있는 기상청의 온도 정보를 이용해서 서울 전체의 온도 지도를 그려보자. 

 

1. 기상청에서 온도 정보를 얻어온다. 

 여기서 귀찮은 수작업이 있는데 오른쪽의 위치를 위도, 경도로 바꿔줘야 한다. 구글 어스로 들어가 위치 정보를 치면 위도, 경도를 구할 수 있다. 

 이렇게 각 위치에 대한 위도 경도를 얻은 후엔 이 정보를 매트랩으로 넘겨준다. 사실 정보라고 해봤자 위도, 경도, 온도 뿐이다. 

데이터 파일

temperature_data_latitude_longitude.mat

2. 이 점들을 데이터로 해서 Gaussian Process Regression을 하자. 

main.m

%% 초기화

clc;

clear all;

close all;

%% 데이터를 불러온다 (서울시의 기상청의 위도 경도와 온도를 담고 있는 데이터)

load 'temperature_data_latitude_longitude.mat';

nr_data = size(z, 1);

data_temp = z(:,1);  

data_lat = z(:,3);  

data_long = z(:,2);

%% 데이터를 약간 가공한다.

% 위치 데이터를 정리한다. 

data_x = (data_lat)*1;

data_y = (data_long)*1;

x_min = min(data_x); x_max = max(data_x);

y_min = min(data_y); y_max = max(data_y);

x_interval = x_max - x_min; y_interval = y_max - y_min;

% 10%의 마진을 주고 x,y의 lim을 구한다.

x_lim = [x_min - x_interval/10 x_max + x_interval/10];

y_lim = [y_min - y_interval/10 y_max + y_interval/10];

% Gaussian Process를 위해서 temperature의 mean을 0으로 한다.

% Gaussian Process는 일반적으로 zero-mean을 가정한다. 

mean_temp = mean(data_temp);

data_temp = data_temp - mean_temp;

%% 데이터를 stem3 plot으로 그려본다.

plot_stem3 = 0; % 그리지 말자. 

if plot_stem3

    stem3(data_x , data_y, data_temp);

    hold on

    plot(data_x, data_y,'ro','MarkerFaceColor',[1 1 1],'MarkerSize',10,'Marker','o',...

        'LineWidth',2);

    hold off

    xlim(x_lim);

    ylim(y_lim);

end

%% GP 로 regresion할 데이터들을 초기화한다.

% 1. GP를 할 위치들이다.

% 우리가 알고 있는 data들을 기반으로 이 새로운 위치에서 온도를 predict한다. 이 때 Gaussian Process Regression을 사용한다. 

% GPR은 우리가 알고 있는 데이터와 새로운 위치가 Gaussian distribution을 갖는다고 가정을 하고 해당 위치에서 온도를 예측. 

res_gp = (min(x_interval, y_interval)/20); % GP 할 때 resolution을 의미한다. (위치 사이의 간격)

x_gp = x_lim(1):res_gp:x_lim(2);

y_gp = y_lim(1):res_gp:y_lim(2);

[x_cord y_cord] = meshgrid(x_gp, y_gp); % 위치를 2차원 grid로 만든다. 

X = [reshape(x_cord, [], 1), reshape(y_cord, [], 1)];

% 2. 데이터를 설정한다. (모두 column vector이다.)

Xbar = reshape([data_x';data_y'], [], 1);

Ybar = data_temp;

%% Gaussian Process Regression을 수행한다. 

sigma2f = 0.3;          % GP kernel의 gain을 의미한다. 

sigma2x = .005;        % GP kernel의 var를 의미한다. 

[zhatxstar sigma2xstar] = GPR_SE_kernel(X, Xbar, Ybar, sigma2f, sigma2x);

%% GPR 결과를 plot한다. 

zhatxstar = zhatxstar + mean_temp; % GPR에 원래 mean을 더해준다. (앞에서 zero mean으로 바꾸기 위해서 평균을 뺐었다.)

option.main_title = 'Gaussian Process Regression on temperature in Seoul @ 2012.11.08';

option.mean_title = 'Temperature mean in Seoul';

option.var_title = 'GPR variance in Seoul';

option.xlabel = 'Longitude';

option.ylabel = 'Latitude';

option.data_x = data_x;

option.data_y = data_y;

mean_gp = reshape(zhatxstar, size(y_gp,2), size(x_gp,2));

var_gp = reshape(diag(sigma2xstar),size(y_gp,2),size(x_gp,2));

figGP = plot_mean_variance(x_gp, y_gp, mean_gp, var_gp, option);

%% 구글 어스에 그려보자. demo_ge_seoul.kml를 실행시켜면 된다. 

addpath('C:\Users\CPSLAB\Dropbox\Research\CPSLab\Implementation\ETC\googleearth');

output = ge_imagesc(x_gp, y_gp, reshape(zhatxstar, size(y_gp,2), size(x_gp,2)));

output2 = ge_colorbar(x_lim(2), y_lim(2), reshape(zhatxstar, size(y_gp,2), size(x_gp,2)));

kmlFileName = 'demo_ge_seoul.kml';

ge_output(kmlFileName, [output2 output], 'name', kmlFileName);

%% 변수 정리 (단위는 위도, 경도, 온도 이다.)

% X: GPR을 수행할 x,y들의 집합 

% Xbar: 측정한 위치 (x, y)

% Ybar: 측정한 값 

% x_gp: x값 만의 벡터 (imagesc에서 사용된다)

% y_gp: y값 만의 벡터 

% mean_gp: 2차원 행렬에 gp로 regression한 mean value들이 들어가있다. 

%           이 값이 이제 reference가 될 것이다. 

% var_gp: 2차원 행렬에 gp로 regression한 varaince value들이 들어있따. 

 

%% 

GPR_SE_kernel.m

function [zhatxstar sigma2xstar] = GPR_SE_kernel(xstar,X,Y,sigma2f,sigma2x)

X = reshape(X',[],1);    

X = reshape(X,2,[])';

sigma2w = 0.3;

n = size(X,1);

KxstarX = GP_K(xstar,X,sigma2f,sigma2x);

KXX = GP_K(X, X, sigma2f, sigma2x);

Kxstarxstar = GP_K(xstar,xstar,sigma2f,sigma2x);

Lambda = (KXX+sigma2w*eye(n))^(-1);

zhatxstar = KxstarX*Lambda*Y;

sigma2xstar = Kxstarxstar - KxstarX * Lambda * KxstarX';

GP_K.m

function K = GP_K(xstar,X,sigma2f,sigma2x)

% Gaussian process 의 Kernel Function이다. 

% 일반적인 SE(Squared exponential) Kernel이다. 

m = size(xstar,1);

c1 = -2*sigma2x;

if(nargin==4)

    n = size(X,1);

    K = zeros(m,n);

    for j=1:n

        for i=1:m

            K(i,j) = exp(norm(xstar(i,:)-X(j,:))^2/c1);

        end

    end

else

    K = eye(m,m);

    for j=1:m-1

        for i=j+1:m

            tmp = exp(norm(xstar(i,:)-xstar(j,:))^2/c1);

            K(i,j) = tmp;

            K(j,i) = tmp;

        end

    end

end    

K = sigma2f* K;

plot_mean_variance.m

function fig = plot_mean_variance(x_gp, y_gp, mean_gp, var_gp, option)

fig = figure('Name',option.main_title...

    ,'Position',[100 100 500 700],'NumberTitle','off'); % figure의 이름과 화면에서 위치, 크기를 설정한다. figure 번호도 없엔다. 

axes('Parent',fig,'Position',[0.12 0.55 0.85 0.39]); 

imagesc(x_gp ,y_gp, mean_gp,'DisplayName','Naive Conditional Mean') % imagesc함수를 이용해서 mean_gp 배열 속에 있는 값들을 그림

colorbar; % 오른족에 각 색이 어떤 값을 나타내는지 그린다. 

title(option.mean_title, 'FontSize', 13); %제목과 글자 크기 

xlabel(option.xlabel);ylabel(option.ylabel); % x,y의 label을 그린다. 

axes('Parent',fig,'Position',[0.12 0.05 0.85 0.39]); % 밑의 그림을 그릴 장소를 정해준다.  

imagesc(x_gp ,y_gp, var_gp,'DisplayName','Naive Conditional Variance') % 마찬가지로 var_gp에 있는 값을 그린다. 

hold on % 현재 그림 위애다가 덮어서 그리겠다는 의미. 

plot(option.data_x, option.data_y,'Color',[1 0 1],'MarkerSize',10,'Marker','x'...

    ,'LineWidth',3,'LineStyle','none','DisplayName','Real Position'); % 하단의 그림 위에다가 option.data_x, option.data_y를 x표시로 그린다.

colorbar; 

title(option.var_title, 'FontSize', 13); % 제목과 글자 크기 

xlabel(option.xlabel);ylabel(option.ylabel);

hold off

이렇게 할 경우 다음과 같은 Gaussian Process 결과를 얻을 수 있다. 이 떄 x, y값은 위도와 경도이다. 

3. 이제 이 정보를 구글 어스 위에 그려보자. 

 구글어스 매트랩 툴박스 (설치하고 적절히 경로 추가)

: addpath('C:\Users\CPSLAB\Dropbox\Research\CPSLab\Implementation\ETC\googleearth');

googleearth_matlab.zip

위의 main을 실행시키면 demo_ge_seoul.kml 파일이 형성된다. 이 파일을 실행하면 다음과 같이 나온다. 

 

서울의 지도 위에 온도 지도가 잘 올라간 것을 알 수 있다.