Skip to content

Instantly share code, notes, and snippets.

@EvangelosAlevizos
Last active February 19, 2020 16:49
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save EvangelosAlevizos/eb164b905f713a215b960e61b91b198b to your computer and use it in GitHub Desktop.
Save EvangelosAlevizos/eb164b905f713a215b960e61b91b198b to your computer and use it in GitHub Desktop.
Here are two matlab scripts that can be used for producing WORLD or KML files for automate georeferencing of images acquired from a UAV.
UAV2WORLD
Input data:
(spatial dimensions expressed in meters)
1) pw, image width (pixels).Right click on the image file click properties select the details tab. insert value within script.
2) ph, image height (pixels).Right click on the image file click properties select the details tab. insert value within script.
3) sw, sensor width (meters).Can be found from the camera specifications. insert value within script.
4) sh, sensor height (meters).Can be found from the camera specifications.insert within script
5) fl, focal length of the camera (meters).Right click on the image file click properties select the details tab. insert value within script.
6) E, easting (meters): X coordinate of the image central pixel
7) N, northing (meters): Y coordinate of the image central pixel
8) alt, altitude (meters above ground)
9) head, heading (positive degrees 0-360)
10) imgs, image file names without the filetype extension:
11) world files extension names should be replaced accordingly within script
Filetype world file extension
GIF .gfw
JPEG .jgw
JPEG 2000 .j2w
PNG .pgw
TIFF .tfw
BMP .bpw or bmpw
-------------------------------------------------------------------------------
UAV2KML
Input data:
(spatial dimensions expressed in decimal degrees)
decimal degrees convertion table: https://en.wikipedia.org/wiki/Decimal_degrees
1) pw, image width (pixels). Right click on the image file click properties select the details tab. insert value within script after converting to decimal degrees.
2) ph, image height (pixels). Right click on the image file click properties select the details tab. insert value within scriptafter converting to decimal degrees.
3) sw, sensor width (decimal degrees).Can be found from the camera specifications. insert value within script after converting to decimal degrees.
4) sh, sensor height (decimal degrees).Can be found from the camera specifications. insert value within script after converting to decimal degrees.
5) fl, focal length of the camera (decimal degrees). Right click on the image file click properties select the details tab. insert value within script after converting to decimal degrees.
6) lon, longitude (decimal degrees): X coordinate of the image central pixel
7) lat, latitude (decimal degrees): Y coordinate of the image central pixel
8) alt, altitude (meters above ground converted to decimal degrees)
9) head, heading (positive degrees 0-360)
10)imgs, image file names with the filetype extension included(e.g.: image1.jpg)
IMPORTANT: After running the script(s) copy and paste the produced files (world or KML) to the directory
where the images are stored. Then, open the georeferenced images either by clicking the KML files
directly (for visualization in Google Earth) or by importing the original images into a GIS program.
%%UAV2KML%%
%% This script outputs one KML file for each image acquired with a UAV for visualization on Google Earth%%
%% Please look at the readme file for detailed instructions about how to use this script%%
%%Evangelos Alevizos, Oct2018
pw=1280; %Image width (pixels)
ph=960; %Image height(pixels)
fl=0.000000006 ; % insert camera focal length from EXIF file and convert in(decimal degrees: dd)
sw=0.000000006332 ; % insert sensor width (dd)
sh=0.000000004725 ; % insert sensor height (dd)
%variables with corresponding matrix names
%Import tables in Matlab--> Image_names: imgs ,longitude(dd): lon ,latitude(dd): lat ,altitude(dd):alt, heading(degrees, 0-360):
%head
for k=1:length(E);
a1=fl/alt(k) ; %scale, dimensionless
a2=1/a1 ; %scale factor
b1=sw/pw ; %dd per pixel on sensor's width
c1=sh/ph ; %dd per pixel on sensor's height
b2=((pw/2)*b1)*a2 ; %offset from central pixel in width direction in dd
c2=((ph/2)*c1)*a2 ; %offset from central pixel in height direction in dd
N1=lat(k)+c2;
S1=lat(k)-c2;
E1=lon(k)+b2;
W1=lon(k)-b2;
%%------KML orientation params----------
Q=cell(17,3);
Q{1,1}='<?xml version="1.0" encoding="UTF-8"?>';
Q{1,2}='';
Q{1,3}='';
Q{2,1}='<kml xmlns="http://www.opengis.net/kml/2.2"> ';
Q{2,2}='';
Q{2,3}='';
Q{3,1}=' <GroundOverlay>';
Q{3,2}='';
Q{3,3}='';
Q{4,1}=' <name>';
Q{4,2}=imgs{k};
Q{4,3}='</name>';
Q{5,1}=' <description></description>';
Q{5,2}='';
Q{5,3}='';
Q{6,1}=' <Icon>';
Q{6,2}='';
Q{6,3}='';
Q{7,1}=' <href>';
Q{7,2}=strcat(imgs{k},'.jpg');
Q{7,3}='</href>';
Q{8,1}=' </Icon>';
Q{8,2}='';
Q{8,3}='';
Q{9,1}=' <LatLonBox>';
Q{9,2}='';
Q{9,3}='';
% NORTH: Specifies the latitude of the north edge of the bounding box, in decimal degrees from 0 to ±90
Q{10,1}=' <north>';
Q{10,2}=num2str(N1,8);
Q{10,3}='</north>';
% SOUTH: Specifies the latitude of the south edge of the bounding box, in decimal degrees from 0 to ±90
Q{11,1}=' <south>';
Q{11,2}=num2str(S1,8);
Q{11,3}='</south>';
% EAST: Specifies the latitude of the east edge of the bounding box, in decimal degrees from 0 to ±180
Q{12,1}=' <east>';
Q{12,2}=num2str(E1,8);
Q{12,3}='</east>';
% WEST: Specifies the latitude of the west edge of the bounding box, in decimal degrees from 0 to ±180
Q{13,1}=' <west>';
Q{13,2}=num2str(W1,8);
Q{13,3}='</west>';
% ROTATION: Specifies a rotation of the overlay about its center, in degrees. Values can be ±180. The default is 0 (north). Rotations are specified in a counterclockwise direction
Q{14,1}=' <rotation>';
Q{14,2}=num2str(head(k));
Q{14,3}='</rotation>';
%
Q{15,1}=' </LatLonBox>';
Q{15,2}='';
Q{15,3}='';
Q{16,1}=' </GroundOverlay>';
Q{16,2}='';
Q{16,3}='';
Q{17,1}='</kml>';
Q{17,2}='';
Q{17,3}='';
%------------------------------------
%writetable(Q,'testKML.txt') ;
%ML=cell2table(Q);
filename = [imgs{k} '.kml'];
fid = fopen(filename,'w');
[rows,cols]=size(Q);
for r=1:rows
fprintf(fid,'%s\t %s\t %s\n',Q{r,:});
%fclose(fid);
end
end
% for k=1:length(E1);
% writetable(Q,'testKML.txt') ;
% % textFilename = ['file' num2str(k) '.txt'];
% % fid = fopen(textFilename, 'w');
% % fprintf(fid,'%d\n %d\n %d\n %d\n %7.2f\n %7.2f\n %s',w1,st);
% % fclose(fid);
% end
%%UAV2WORLD%%
%% This script outputs one world file for each image acquired with a UAV for automate georeferencing%%
%% Please look at the readme file for detailed instructions about how to use this script%%
%%Evangelos Alevizos, Oct2018
pw=1280; %insert image width from EXIF file(pixels)
ph=960; %insert image heightfrom EXIF file(pixels)
fl=0.006 ; %insert focal length from EXIF file(meters)
% sw=0.006332 ; %insert sensor width (meters)
% sh=0.004725 ; %insert sensor height (meters)
%variables
%E,N,alt
for z=1:length(E);
a1=fl/alt(z) ; %scale, dimensionless
a2=1/a1 ; %scale factor
b1=sw/pw ; %meters per pixel on sensor's width
c1=sh/ph ; %meters per pixel on sensor's height
b2=((pw/2)*b1)*a2 ; %offset from central pixel in width direction in meters
c2=((ph/2)*c1)*a2 ; %offset from central pixel in height direction in meters
%
% world file structure (row-wise): cellsize_x,rotation param_x,rotation param_y,-cellsize_y,upperleft_E,upperleft_N,
%%Import tables in Matlab--> Image_filenames: imgs ,Easting(meters): E ,Northing(meters): N ,altitude(meters):alt, heading(degrees, 0-360):
%head
czx=b1*a2 ; %cellsize in width direction in meters
czy=(c1*a2) ; %cellsize in width direction in meters
E1=E(z)+(b2) ; %upper left pixel's X-coordinate in meters
N1=N(z)-(c2) ; %upper left pixel's y-coordinate in meters
ca = cos(degtorad(head(z))) * czx;
cb = -(sin(degtorad(head(z))) * czy);
cc = -(sin(degtorad(head(z))) * czx);
cd = -(cos(degtorad(head(z))) * czy);
% zxn= sqrt((ca^2)+(cb^2));
% zyn= sqrt((cc^2)+(cd^2));
%for k=1:length(E);
w1=[ca;cb;cc;cd;E1;N1] ; %world file data
textFilename = [imgs{z} '.bpw'];
fid = fopen(textFilename, 'w');
fprintf(fid,'%d\n %d\n %d\n %d\n %7.2f\n %7.2f\n',w1);
fclose(fid);
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment