How can I convert longitude and latitude in x and y coordinates on a Mollweide map?
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
I have to display the coordinates attached to the question on a Mollweide map. To do so I've tried implementing the Newton iterative method shown in the wikipedia page "Mollweide projection", like so:
R = 6.378145*(10^6);
axis([-2*R*sqrt(2) 2*R*sqrt(2) -R*sqrt(2) R*sqrt(2)]);
told = latitude(1);
t = told;
teta = nan(length(latitude), 1);
x = nan(length(longitude), 1);
y = nan(length(latitude), 1);
% Iterative cycle
for i = 1 : length(latitude)
while (abs(2*t + sin(2*t) - pi*sin(latitude(i))) >= 1e-6)
tnew = told - (2*told+sin(2*told)-pi*sin(latitude(i)))/(2+2*cos(2*told));
told = tnew;
t = tnew;
end
teta(i) = t;
x(i) = (R*2*sqrt(2)/pi)*longitude(i)*cos(teta(i));
y(i) = R*sqrt(2)*sin(teta(i));
end
x = x*pi/180;
Now this doesn't seem to work since the attached points, that can be seen by running this simple code:
figure(2);
hold on;
currIdx = 1;
for i = [find(abs(diff(longitude')) > 10), length(longitude)]
idx = currIdx:i;
if length(idx) < 2
style = '.r';
else
style = 'r';
end
plot(longitude(idx), latitude(idx), style);
axis([-180 180 -90 90]);
currIdx = i+1;
end
defenitely don't match with those obtained with the first piece of code. They can be visualized in the same way:
figure(3);
axis([-2*R*sqrt(2) 2*R*sqrt(2) -R*sqrt(2) R*sqrt(2)]);
hold on;
currIdx = 1;
for i = [find(abs(diff(x')) > 1.5e7-1e7), length(y)]
idx = currIdx:i;
if length(idx) < 2
style = '.r';
else
style = 'r';
end
plot(x(idx), y(idx), style);
currIdx = i+1;
end
Why does this happen and how can I get the correct coordinate conversion?
1 commentaire
Christian
le 15 Août 2019
I need to do the same thing. I hope you don't mind me asking before I get into it myself, but did you find any solution up to now?
Réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!