Mapping - latitude-longitude / projection of a circle

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

Answer 1

[Meg Noah]

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 
earthRadius_m = 6371.23e3;

% 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 
beta_deg = rad2deg(distance_m/earthRadius_m);

% 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))']);
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 = 'http://upload.wikimedia.org/wikipedia/commons/thumb/c/cd/Land_ocean_ice_2048.jpg/1024px-Land_ocean_ice_2048.jpg';
cdata = imread(image_file);
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');

[darova]

Finally! Took 9 years! Where have you been all this time?

[Meg Noah]

:-)

I just have a couple of weekends before the next semester starts. Wanted to see what the 'Answer' section is all about, and how well it works.

[devadas kuna]

hi sir...

in the above code "fov_deg = 30; " what it indicates..

and latitude0_degN = 20; % [degN] nadir intersection (center of FOV)

longitude0_degE = -75.2; % [degE] nadir intersection (center of FOV)

is they are geodetic receiver lat,long??

[Meg Noah]

FOV = Field of View. That's the angular extent of either the horizontal or vertical projection into the receiver.

In the solution, the latitude and longitude are for a spherical earth - polar coordinates but with 'colatitude' instead of the typical polar angle from the z-axis. You need the mapping toolbox or to write the other functions for geodetic coordinates. If you need that I can provide you the math for the geodetic conversions to WGS84, but why not just use the matlab routines?

[zeynep ozkayikci]

when i try to zoom the figure, the picture is distorted. How can i fix it ? do you have an idea, thanks

[Meg Noah]

Hi,

It is possible. From the dropdown menu at the top of the figure, select:

View | Camera Toolbar

Click the 5th button (5th from left side), to Move Camera Forward and Backward

[zeynep ozkayikci]

Thank you so much, you make me so happy :)