Mapping - latitude-longitude / projection of a circle

45 views (last 30 days)
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.

Answers (1)

Meg Noah
Meg Noah on 10 Jan 2020
Here's a simple solution
close all
clear all
% select case you want
icase = 2; % 1=landat7 width of image, 2=geosynchronous
switch icase
case 1
% 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)
% approximate with spherical earth
earthRadius_m = 6371.23e3;
% from the law of cosines
distance_m = (altitudeSatellite_m)*acos(sind(latitude0_degN)*sind(latitude0_degN) + ...
% 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
beta_deg = rad2deg(distance_m/earthRadius_m);
% distance + heading formulae
latitude_degN = asind(sind(latitude0_degN).*cosd(beta_deg) + ...
longitude_degE = longitude0_degE + ...
atan2d(sind(angle_deg).*sind(beta_deg).*cosd(latitude0_degN), ...
% if you have aerospace toolbox - can do ellipsoidal earth
% posECEF_m = lla2ecef([latitude_degN',longitude_degE',zeros(size(latitude_degN))']);
x = earthRadius_m*cosd(latitude_degN).*cosd(longitude_degE);
y = earthRadius_m*cosd(latitude_degN).*sind(longitude_degE);
z = earthRadius_m*sind(latitude_degN);
% draw 3D earth
npanels = 60;
alpha = 0.5; % globe transparency level, 1 = opaque, through 0 = invisible
image_file = '';
cdata = imread(image_file);
hold on; set(gca, 'NextPlot','add', 'Visible','off');
axis equal; axis auto;
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

Sign in to comment.


Find more on Coordinate Systems in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!