Azzera filtri
Azzera filtri

Mapping - latitude-longitude / projection of a circle

19 visualizzazioni (ultimi 30 giorni)
Geant Bepi
Geant Bepi il 21 Feb 2011
Commentato: zeynep ozkayikci il 11 Ago 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

Risposte (1)

Meg Noah
Meg Noah il 10 Gen 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
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');
  7 Commenti
Meg Noah
Meg Noah il 11 Ago 2022
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
zeynep ozkayikci il 11 Ago 2022
Thank you so much, you make me so happy :)

Accedi per commentare.

Categorie

Scopri di più su Reference Applications in Help Center e File Exchange

Prodotti

Community Treasure Hunt

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

Start Hunting!

Translated by