Thread Subject: lats and lons

Hi, I posted earlier. I thought I would repost with more
clearness as I still have not solved the issue:)
So I have this matrix subdat_sat_interp of size 53x71.
I plot it in the following way:

h=pcolor(subdat_sat_interp(:,:)'); hold on; axis xy;

I also have two longitudes and latitudes arrays (lon_loc and
lat_loc) of size 53x71, which I plot on the same graph
using countour function as:

[c,h]=contour(lon_loc');
clabel(c,h);
[c,h]=contour(lat_loc');
clabel(c,h)

Note, lon_loc ranges from [0.0565, 359.9424]
And lat_loc ranges from [5.1660, 90.0].

Now I have two vectors (la and fi) of size 1x137 that have
latitudes and longitudes respectively ranging from [57.2,
78.3] and [5.3, 358.8] respectively. The goal for me now is
to plot all of these 137 la and fi on the same graph as
previously (see above) as points. So I want to plot:
plot (x,y,'or'). where x and y should be of size 1x137
ranging from [1,71] and [1,53] respectively. So I need some
function that would transfer my la and fi vectors into
respective x and y. Any ideas? THANKS!!!!!

 > I also have two longitudes and latitudes arrays (lon_loc and
> lat_loc) of size 53x71, which I plot on the same graph
> using countour function as:
>
> [c,h]=contour(lon_loc');
> clabel(c,h);
> [c,h]=contour(lat_loc');
> clabel(c,h)
>
> Note, lon_loc ranges from [0.0565, 359.9424]
> And lat_loc ranges from [5.1660, 90.0].
>
> Now I have two vectors (la and fi) of size 1x137 that have
> latitudes and longitudes respectively ranging from [57.2,
> 78.3] and [5.3, 358.8] respectively. The goal for me now is
> to plot all of these 137 la and fi on the same graph as
> previously (see above) as points. So I want to plot:
> plot (x,y,'or'). where x and y should be of size 1x137
> ranging from [1,71] and [1,53] respectively. So I need some
> function that would transfer my la and fi vectors into
> respective x and y. Any ideas? THANKS!!!!!

I'm assuming that you have something similar to projected
data, but without knowing the specific details of the
projection. You should be able to "project" your data
points using griddata with your lat and lon grids. The
example below demonstrates this (I used some mapping toolbox
functions to generate my fake data, but the toolbox isn't
needed for the projection part).

% Create some lat/lon gridded data (your lon_loc and lat_loc)

axesm('robinson');
x = linspace(-1.8, -.8, 53);
y = linspace(.3, 1.1, 71);
[x,y] = meshgrid(x,y);
[lat, lon] = minvtran(x, y);
close;

[x,y] = meshgrid(1:53, 1:71);

% Points to plot (your la and fi)

S = load('coast.mat');
latvec = S.lat;
lonvec = S.long;

% Figure out corresponding x and y coords

xvec = griddata(lat, lon, x, latvec, lonvec)
yvec = griddata(lat, lon, y, latvec, lonvec);

% Plot lat and lon contours with "projected" points

figure;
[c1,h1]=contour(x,y, lon);
clabel(c1,h1);
hold on;
[c2,h2]=contour(x,y,lat);
clabel(c2,h2);
plot(xvec,yvec, '-b.');

Kelly, thank you soooooooooo much. This is exactly what I
have been looking for. THANKS!!!


"Kelly Kearney" <kakearney@nospamgmail.com> wrote in message
<g827e9$jik$1@fred.mathworks.com>...
> > I also have two longitudes and latitudes arrays
(lon_loc and
> > lat_loc) of size 53x71, which I plot on the same graph
> > using countour function as:
> >
> > [c,h]=contour(lon_loc');
> > clabel(c,h);
> > [c,h]=contour(lat_loc');
> > clabel(c,h)
> >
> > Note, lon_loc ranges from [0.0565, 359.9424]
> > And lat_loc ranges from [5.1660, 90.0].
> >
> > Now I have two vectors (la and fi) of size 1x137 that have
> > latitudes and longitudes respectively ranging from [57.2,
> > 78.3] and [5.3, 358.8] respectively. The goal for me now is
> > to plot all of these 137 la and fi on the same graph as
> > previously (see above) as points. So I want to plot:
> > plot (x,y,'or'). where x and y should be of size 1x137
> > ranging from [1,71] and [1,53] respectively. So I need some
> > function that would transfer my la and fi vectors into
> > respective x and y. Any ideas? THANKS!!!!!
>
> I'm assuming that you have something similar to projected
> data, but without knowing the specific details of the
> projection. You should be able to "project" your data
> points using griddata with your lat and lon grids. The
> example below demonstrates this (I used some mapping toolbox
> functions to generate my fake data, but the toolbox isn't
> needed for the projection part).
>
> % Create some lat/lon gridded data (your lon_loc and lat_loc)
>
> axesm('robinson');
> x = linspace(-1.8, -.8, 53);
> y = linspace(.3, 1.1, 71);
> [x,y] = meshgrid(x,y);
> [lat, lon] = minvtran(x, y);
> close;
>
> [x,y] = meshgrid(1:53, 1:71);
>
> % Points to plot (your la and fi)
>
> S = load('coast.mat');
> latvec = S.lat;
> lonvec = S.long;
>
> % Figure out corresponding x and y coords
>
> xvec = griddata(lat, lon, x, latvec, lonvec)
> yvec = griddata(lat, lon, y, latvec, lonvec);
>
> % Plot lat and lon contours with "projected" points
>
> figure;
> [c1,h1]=contour(x,y, lon);
> clabel(c1,h1);
> hold on;
> [c2,h2]=contour(x,y,lat);
> clabel(c2,h2);
> plot(xvec,yvec, '-b.');
>

"jenya polyakova" <jenya56@yahoo.com> wrote in message
<g829nd$d18$1@fred.mathworks.com>...
> Kelly, thank you soooooooooo much. This is exactly what I
> have been looking for. THANKS!!!

Using gridddata() has *got* to be a bad solution for this
problem when the real latitudes range up to 78.3 degrees.
At 78.3 degrees, a degree of longitude translates into
only roughly 14 miles instead of the 69-some-odd that
it would be at sea-level at the equator. With the ranges
of latitudes and longitudes you are using, you have
significant non-linear projection distortion, which the
griddata() solution is *not* taking into account!

The longitudes you mention for the containing grid
are very nearly spanning the entire world, so unless you
have proof otherwise, you should be assuming that some
of the paths you want to plot might wrap around the edges
of the grid rather than staying within the grid. And
with that big of a range of lats and longs you really
should be plotting great circles instead of approximating
by straight lines.

could you please suggest how would I then (without using
griddata) approach my problem? I tried several other ways
but with no success. Thanks.
           
"Walter Roberson" <roberson@ibd.nrc-cnrc.gc.ca> wrote in
message <g82bqt$7vi$1@fred.mathworks.com>...
> "jenya polyakova" <jenya56@yahoo.com> wrote in message
> <g829nd$d18$1@fred.mathworks.com>...
> > Kelly, thank you soooooooooo much. This is exactly what I
> > have been looking for. THANKS!!!
>
> Using gridddata() has *got* to be a bad solution for this
> problem when the real latitudes range up to 78.3 degrees.
> At 78.3 degrees, a degree of longitude translates into
> only roughly 14 miles instead of the 69-some-odd that
> it would be at sea-level at the equator. With the ranges
> of latitudes and longitudes you are using, you have
> significant non-linear projection distortion, which the
> griddata() solution is *not* taking into account!
>
> The longitudes you mention for the containing grid
> are very nearly spanning the entire world, so unless you
> have proof otherwise, you should be assuming that some
> of the paths you want to plot might wrap around the edges
> of the grid rather than staying within the grid. And
> with that big of a range of lats and longs you really
> should be plotting great circles instead of approximating
> by straight lines.

Walter wrote:

> Using gridddata() has *got* to be a bad solution for this
> problem when the real latitudes range up to 78.3 degrees.
> At 78.3 degrees, a degree of longitude translates into
> only roughly 14 miles instead of the 69-some-odd that
> it would be at sea-level at the equator. With the ranges
> of latitudes and longitudes you are using, you have
> significant non-linear projection distortion, which the
> griddata() solution is *not* taking into account!


Well, I did make several assumptions about her data. First, the
latitude range she originally mentioned did not span the entire
northern hemisphere, only about 55 N to 80 N, which reduces the range
of distortion a bit. And though she didn't mention one way or the
other, the fact that she begins with a rectangular grid for which a
latitude and longitude value are provided for each grid cell made me
assume she was dealing with a cylindrical projection that projects to
a rectangular map, such as a miller or mercator projection. I also
assumed, based on her range of longitudes from 5 to 358 that none of
the objects she wanted to plot crossed the prime meridian.

If this is the case, the griddata method performs very well, and does
replicate the expected exaggeration at the poles. It will likely fail
if she trying to plot very small objects (e.g. a shape that falls
entirely within one of her original grid cells).

Of course, if she's dealing with a projection that's not close to
rectangular at high latitudes (a sinusoidal projection, for example),
then it will probably fail miserably. :-)

Maybe the M_MAP toolbox by Rich Pawlowicz could help.

http://www.eos.ubc.ca/~rich/#M_Map

Carlos

jenya polyakova wrote:
> could you please suggest how would I then (without using
> griddata) approach my problem? I tried several other ways
> but with no success. Thanks.

The reason you have had no success so far is that you have
not clearly defined what you are trying to do.

As I wrote in one of your other threads:

>>>Latitudes and longitudes are not enough information to convert
>>>a planar grid. Latitudes and longitudes give you two coordinates
>>>in a spherical coordinate system, but you need the third coordinate
>>>(radius) in order to be able to use the information. Without
>>>the third coordinate, and knowledge of the projection system in use
>>>(e.g., Mercator, Robinson, Lambert) you cannot project onto a
>>>planar grid. Assuming a constant radius is *not* a good approximation
>>if a significant angular latitude range is being used.

And sure enough, you -are- using a significant angular latitude range.
You are not just trying to create a map of a few square miles, where
the spherical distortion can be ignored.

You have not clearly defined your coordinate system, whether it is
geographic or geocentric, and you have not taken into account the
oblateness of the Earth. On the other hand, according to the constants at
http://topex.ucsd.edu/geodynamics/14gravity1_2.pdf
and a bit of calculation, the maximum geographic vs geocentric
latitude distance works out to be roughly 12 seconds, which is
not significant on the coarse grid you are using. Likewise, the
oblateness is roughly one part in 300; as your grid is only 71 units
for latitude, I suspect you won't care if a point very near a grid
boundary trickles over to the next element.

Mostly, though, you have not defined your projection system and
any accompanying altitude data (or the reference point(s) that
the altitudes are defined relative to.) Saying that you want "a grid"
is not sufficient. *Are* you assuming a cylindrical projection? And do
you care about the significant distortions that is going to introduce
towards the north of your graph? Or to phrase it a different way:
this graph that you are producing, what are it's important characteristics?
Do the horizontal and vertical need to be on the same scale so that you
can do measurements in arbitrary directions? Do the horizontal and
vertical each need to be on linear scales so that a distance measured
in a fixed direction near the top of the graph would be comparable to
a distance measured in the same direction near the bottom of the graph?
Is it important to be able to measure angles off of the graph?

Or is this whole thing just about creating a vague web-page icon of
the Northern Hemisphere, where correctness in any mathematical sense
is not important as long as it looks more or less okay artistically?

--
Q = quotation(rand);
if isempty(Q); error('Quotation server filesystem problems')
else sprintf('%s',Q), end

Ok, if you just want a quick and nasty view of how a series
of lat longs look relative to each other avoiding the issue
of projections, then assume a perfectly spherical earth and
use sph2cart to convert to a set of points:

[x,y,z] = sphere(50);

% plot sphere
surf(x,y,z)

% create some lat/longs
lat = linspace(-32.215, -32.15, 100);
lon = linspace(115.15, 115.22, 100);

[px py pz] = sph2cart(lon, lat, radius);

axis equal
hold on
plot3(px, py, pz, 'k.')

Rod

--
http://iheartmatlab.blogspot.com

Rodney Thomson wrote:
> Ok, if you just want a quick and nasty view of how a series
> of lat longs look relative to each other avoiding the issue
> of projections, then assume a perfectly spherical earth and
> use sph2cart to convert to a set of points:

> [x,y,z] = sphere(50);

Interestingly, if this is done, then the choice of radius is
not important: the radius just becomes a uniform scale factor
on the graph, and using different radii would just correspond
to using different units. 6378137 meters, 3963.3 miles,
35.5 kizlodes... doesn't matter unless you want to take
a measurement from the graph.

--
Q = quotation(rand);
if isempty(Q); error('Quotation server filesystem problems')
else sprintf('%s',Q), end