% From latitude longitude geographical coordinate system to local planar 
% coordinate system involves translation, rotation, and scaling of 
% coordinate system. 
%Input lat,lon of the points; 
% (lat0,lon0):	the latitude longitude for the origin point of the local planar coordinate system 
% azimuth: the azimuth angle of x-axis of the local coordinate system, 
% y-axis is clockwise from x-axis. 
% 
% output x=xy(2); y=xy(1); 
% (x,y): the points location of local planar coordinate system. unit in meter 
% 
% Notice: Not suitable for calculing large scale points. For points within 1 km, 
% the accuracy depends on the azimuth angle.

% example: Set the local origin point (0,0) at (55.709264 N, 13.201449 E), a building corner 
% the azimuth angle of x-axis is 104.5 deg 
% the the point 55.709224,13.202661) of the geographic coordinate 
% system has local coordinates at 
% xy = latlon2xy(55.70922447,13.20266068,55.709264,13.201449,104.5); 
% 
% See also xy2latlon 
% Hongxiao Jin Nov 29 2014 
% 
function xy = latlon2xy(lat,lon,lat0,lon0,azimuth) 
xy = []; 
if size(lat,1) >size(lat,2), lat=lat'; end 
if size(lon,1) >size(lon,2), lon=lon'; end 
if size(lat,1) ~= 1 
msgbox('latitude should be a vector!');return; 
end 
if size(lon,1) ~= 1 
msgbox('longitude should be a vector!');return; 
end 
if size(lat,2) ~= size(lon,2) 
msgbox('Latitude Longitude should be vectors of equal length !');return; 
end

R = 6371007.181; % earth radius Spheric model, same as Sinosoidal projection system 
PI=3.141592653589793238; 
theta=azimuth-90; 
xy = [cos(PI*theta/180) sin(PI*theta/180); -sin(PI*theta/180) cos(PI*theta/180)]*[1 0;0 cos(PI*lat0/180)]*([lat;lon] - [lat0 lon0]'*ones(size(lat))) *PI*R/180;