# Mapping - latitude-longitude / projection of a circle

64 views (last 30 days)
Geant Bepi on 21 Feb 2011
Commented: zeynep ozkayikci on 11 Aug 2022
Ok! I'll try to break the problem in to simpler parts.
There's a circle above a sphere ( for e.g., a view-cone of a satellite orbiting the Earth). Plot the area covered by the circle (so, essentially a small-circle with a radius "r" (scircle1)) on a latitude-longitude map. Consider an orbit around a perfect sphere.
We know the location of the satellite (i.e., latitude (theta), longitude(phi), altitude) and the radius of the circle (alfa). Throughout the calculation longitude and altitude remain constant.
The position and attitude of the five windows (above mentioned view-cone) of the satellite looking at different directions (one is pointed upwards and the other four are perpendicular to the one pointed upwards and perpendicular to each other) is calculated from the following rotation matrix;
R1=[cos(theta) 0 sin(theta); 0 1 0; -sin(theta) 0 cos(theta)]
R2=[cos(pi/4) sin(pi/4) 0; -sin(pi/4) cos(pi/4) 0; 0 0 1]
from R1R2 matrix we get five unit vectors; n1=[sin(atheta1); 0; cos(atheta1)]; n2=[cos(atheta1)/sqrt(2); -1/sqrt(2); -sin(atheta1)/sqrt(2)]; n3=[cos(atheta1)/sqrt(2); 1/sqrt(2); sin(atheta1)/sqrt(2)]; n4=-n2; n5=-n3;
My problem starts here. For n1 I used a different method and drew the area covered by the circle as follows;
Since we know latitude, longitude and radius of the circle,
[atheta1,aphi1]=antipode(theta1,phi1); [theta2,phi2]=scircle1(atheta1,aphi1,alfa); H = plot(phi2,theta2,'--k');
This part of the code does exactly what I want and plot the boundary of the area covered by the first circle (represented by n1), on the latitude-longitude map that I am generating. It is important to note that our final projection should be on this latitude-longitude map.
But I have no clue as to how to repeat this to n2 through n5.
Probably you already noticed that we did not need any help from the rotation matrix or n1 unit vector to draw the area covered by the first circle. Since we have five circles I calculated the above five unit vectors (n1 - n5) that represent each of the circles.
How do I move forward from here. How do I get MATLAB to solve n1 so that I can plot the area covered by it.
I hope i explained the problem clearly. Any tip, help is greatly appreciated.
Thanks

Meg Noah on 10 Jan 2020
Here's a simple solution
clc
close all
clear all
% select case you want
icase = 2; % 1=landat7 width of image, 2=geosynchronous
switch icase
case 1
% https://earth.esa.int/web/eoportal/satellite-missions/l/landsat-7
% ifov = 42.5 microradians for a single detector with ~30 m resolution
% swath width is 15deg FOV (~185 km) -> 6160 pixels
fov_deg = 6160*rad2deg(42.5e-6); %
altitudeSatellite_m = 705.0e3; % [m] sun-synchronous polar orbit
% this solution only works for NADIR pointing satellite FOV
% location of center of footprint
latitude0_degN = 42; % [degN] nadir intersection (center of FOV)
longitude0_degE = -72; % [degE] nadir intersection (center of FOV)
case 2
altitudeSatellite_m = 35786e3; % [m] geosynchronous
fov_deg = 30; %
% this solution only works for NADIR pointing satellite FOV
% location of center of footprint
latitude0_degN = 20; % [degN] nadir intersection (center of FOV)
longitude0_degE = -75.2; % [degE] nadir intersection (center of FOV)
end
% approximate with spherical earth
% from the law of cosines
distance_m = (altitudeSatellite_m)*acos(sind(latitude0_degN)*sind(latitude0_degN) + ...
cosd(latitude0_degN)*cosd(latitude0_degN)*cosd(fov_deg));
% sweep out distances from lat0, lon0 at distance_m through bearings 0:360
angle_deg = 0:360;
% earth centered angle swept out by the FOV
% distance + heading formulae
latitude_degN = asind(sind(latitude0_degN).*cosd(beta_deg) + ...
cosd(latitude0_degN).*sind(beta_deg).*cosd(angle_deg));
longitude_degE = longitude0_degE + ...
atan2d(sind(angle_deg).*sind(beta_deg).*cosd(latitude0_degN), ...
cosd(beta_deg)-sind(latitude0_degN).*sind(latitude_degN));
% if you have aerospace toolbox - can do ellipsoidal earth
% posECEF_m = lla2ecef([latitude_degN',longitude_degE',zeros(size(latitude_degN))']);
% draw 3D earth
npanels = 60;
alpha = 0.5; % globe transparency level, 1 = opaque, through 0 = invisible
figure('color','white');
hold on; set(gca, 'NextPlot','add', 'Visible','off');
axis equal; axis auto;
view(0,30);
axis vis3d;
[xe, ye, ze] = ellipsoid(0, 0, 0, earthRadius_m, earthRadius_m, earthRadius_m, npanels);
globe = surf(xe, ye, -ze, 'FaceColor', 'none', 'EdgeColor', 0.5*[1 1 1]);
set(globe, 'facecolor', 'texturemap', 'cdata', cdata, 'facealpha', ...
alpha, 'edgecolor', 'none');
% annotate footprint as a circle
plot3(x,y,z,'r');
zeynep ozkayikci on 11 Aug 2022
Thank you so much, you make me so happy :)