All times are UTC  Page 2 of 3 [ 33 posts ] Go to page Previous  1, 2, 3  Next
 Print view Previous topic | Next topic
Author Message
 Post subject: Re: LRO - WAC DEM Posted: Sun, 06-10-13, 16:28 GMT  Joined: Fri, 03-04-09, 8:21 GMT
Posts: 221
cartrite wrote:
I got this output. value = -6479 Descaled Value = 1734160.5 note* not sure what Descaled Value is?
Possible first line of tab file. The lat/long coords are probably correct. The value still questionable.
Code:
-86.46504 -45.00000 1730.921
If you notice the line in your gdalinfo output:
Quote:
Offset: 1737400, Scale:0.5
1737400 is the nominal reference radius of the Moon in metres, 0.5 is the scale factor to convert the values to meters. (Both of these info's are also in the LBL file.)
The descaled value (I would guess) is the elevation in meters from the center of the Moon (rather than the surface).
ie. by the following formula:
Code:
Descaled value = Offset + ( value * scale factor )
= 1737400 + (-6479 * 0.5)
= 1737400 - 3239.5
= 1734160.5
... in other words, 3239.5 m below the nominal surface.

The lats and longs you report do seem to be in agreement with your comments regarding the corners being outside of (north of) the 87.5 deg latitude line, which is just further confirmation that the raw file coordinates are cartesian in nature.

Top Post subject: Re: LRO - WAC DEM Posted: Sun, 06-10-13, 17:04 GMT  Joined: Thu, 25-10-07, 15:20 GMT
Posts: 992
Location: NE PA, USA
I've been playing around with gdal all morning and haven't come up with a quick answer to creating a lat long file. I did come across this file. http://pds-geosciences.wustl.edu/lunar/ ... _polar.cat

It has formulas to get lat long coordinates.

Quote:
The spherical equations relating map coordinates (x, y) to
planetocentric coordinates (Lat, Lon) [SNYDER1987] are:

North Polar Stereographic
x = 2 * Rp * TAN(Pi / 4 - Lat / 2) * SIN(Lon - LonP)
y = -2 * Rp * TAN(Pi / 4 - Lat / 2) * COS(Lon - LonP)

South Polar Stereographic
x = 2 * Rp * TAN(Pi / 4 + Lat / 2) * SIN(Lon - LonP)
y = 2 * Rp * TAN(Pi / 4 + Lat / 2) * COS(Lon - LonP)

Where LonP is the central longitude, LatP is the latitude of true
scale and is always 90 or -90, and Rp is the polar radius of the Moon.
LonP corresponds to the +y axis, i.e., zero longitude is up in the map.

The spherical inverse formulas for Lat and Lon from X and Y position
in the image array are:

Lat = ARCSIN[COS(C) * SIN(LatP) + y * SIN(C) * COS(LatP) / P]

North Polar Stereographic
Lon = LonP + ARCTAN[x / (-y)]

South Polar Stereographic
Lon = LonP + ARCTAN[x / y]

where:
P = SQRT(x^2 + y^2)
C = 2 * ARCTAN(P / 2 * Rp)

recall:
x = (Sample - S0) * Scale
y = (L0 - Line) * Scale

This could probably be used to write a program or script to generate the coords.

LatP would be -90
LonP would be 0
Rp would be 1737400
etc.
cartrite

Top Post subject: Re: LRO - WAC DEM Posted: Sun, 06-10-13, 19:21 GMT  Joined: Tue, 04-09-07, 21:55 GMT
Posts: 873
Location: N 42.38846 W 83.45456
Quote:
so I'm assuming that all the data contained within the file describes only the circle of radius 2.5 degrees from the south pole, rather than a rectangular swathe of geography centered on the pole, and therefore when we present it as a rectangular image, we're actually stretching it's true geography in the corners.

No
the corners are higher lat of the image
it is just NOT cropped

Quote:
I'm not trying to do a complete spherical moon, but rather the idea is to place a patch of geometry at the polar region. (Same idea as what Selden did with his Mt Palomar telescope, IIRC).

i would use isis , but GDAL should be able to crop it ( in isis it is putting a check in a check box and running it to crop out the extra px and make it a circle )

the curvature from the center to the edge is not going to be much

i am on my first cup of coffee for the day so -- no math

calculate the height of a part of the sphere with a plane bisecting at -87.5 and the center at -90

and use the gmic to .off mesh code

in the very last post , right above this
Code:
MAP_RESOLUTION = 1516.17 <pix/deg>

when i was remapping it to work with a smaller image
at 128 ppd it is a image that is 641x641 px

Top Post subject: Re: LRO - WAC DEM Posted: Tue, 08-10-13, 13:01 GMT  Joined: Fri, 03-04-09, 8:21 GMT
Posts: 221
Progress...

I can now generate a CMOD directly from the IMG file. Still some imperfections and issues to resolve, but it is progress.

Here are some overviews of the generated terrain in place at the South Pole along with my Shackleton Base addon (CLICK for fullsize)....

Attachment:

Attachment:

Attachment:

Last edited by chuft-captain on Tue, 08-10-13, 13:09 GMT, edited 2 times in total.

Top Post subject: Re: LRO - WAC DEM Posted: Tue, 08-10-13, 13:07 GMT  Joined: Fri, 03-04-09, 8:21 GMT
Posts: 221
And some additional views (CLICK for fullsize)....

Shackleton Base up close:
Attachment:

A view from high up inside the dome overlooking buildings far below and the lunar terrain outside:
Attachment:

Shackleton Base and terrain as it appears from the Malapert Mountains communication centre:
Attachment:

Top Post subject: Re: LRO - WAC DEM Posted: Tue, 08-10-13, 16:28 GMT  Joined: Thu, 25-10-07, 15:20 GMT
Posts: 992
Location: NE PA, USA
chuft-captain wrote:
Progress...

I can now generate a CMOD directly from the IMG file. Still some imperfections and issues to resolve, but it is progress.

How did you end up generating the cmod?

cartrite

Top Post subject: Re: LRO - WAC DEM Posted: Wed, 09-10-13, 1:03 GMT  Joined: Fri, 03-04-09, 8:21 GMT
Posts: 221
I extended my original "script" (which generated the colorised image map) to also generate the CMOD.
I made use of some of the formulae you provided above, and I have replicated and customized some of Fridger's logic from his PERL script for the tessellation.
There's a few subtle differences, and it's still an imperfect first cut, so I may yet need to consult with The Doctor on the finer details of his logic. To avoid generating a 1.2 Gigabyte CMOD, I'm utilizing only every 10th pixel to create a 758x758 mesh (which is still 75MB in size). This means that the resolution of the terrain is only 200m between samples rather than the 20m/pix resolution of the image.

Top Post subject: Re: LRO - WAC DEM Posted: Wed, 09-10-13, 1:28 GMT  Joined: Thu, 25-10-07, 15:20 GMT
Posts: 992
Location: NE PA, USA
There is a smaller file in the folder. One at 240 meters a pixel. It extends to -75 lat though. It may be worth it though because the values are probably interpolated. It can always be cropped if it is too big.

I tried using the formulas to write a program that creates a (lat / lon / height) file like the phobos.tab but it doesn't work. It can do the longitudes but the latitudes are way off. Probably something to do with radians and degrees. Trig functions in c++ return all the values in radians. I think they are getting mixed somehow.

cartrite

Top Post subject: Re: LRO - WAC DEM Posted: Wed, 09-10-13, 1:41 GMT Site Admin Joined: Fri, 31-08-07, 7:01 GMT
Posts: 4675
Location: Hamburg, Germany
cartrite wrote:
There is a smaller file in the folder. One at 240 meters a pixel. It extends to -75 lat though. It may be worth it though because the values are probably interpolated. It can always be cropped if it is too big.

I tried using the formulas to write a program that creates a (lat / lon / height) file like the phobos.tab but it doesn't work. It can do the longitudes but the latitudes are way off. Probably something to do with radians and degrees. Trig functions in c++ return all the values in radians. I think they are getting mixed somehow.

cartrite

Converting degrees in radians and vice versa is very simple. I am sure you know that.
If not

So just use angle_radian throughout in trig functions in C++, Perl code. Never mix!

Good luck,
Fridger

Top Post subject: Re: LRO - WAC DEM Posted: Wed, 09-10-13, 2:25 GMT  Joined: Thu, 25-10-07, 15:20 GMT
Posts: 992
Location: NE PA, USA
This was just an exercise that went blewie. I was just trying to see if I could create the map coordinates.
My last attempt looked like this.
Code:
// ARCSIN[COS(C) * SIN(LatP) + y * SIN(C) * COS(LatP) / P]
//                lat = asin (cos(C) * sin(cen_lat) + (y * sin(C) * cos(cen_lat) / P)) * 180.0 / PI;

The top line was the formula from http://pds-geosciences.wustl.edu/lunar/ ... _polar.cat
The bottom line was what I tried. It produced a list of latitudes way out of range.
The closest try was when I tried this. This was still out of range.
Code:
lat = asin (cos(C * 180.0 / PI) * sin(cen_lat * 180.0 / PI) + (y * sin(C * 180.0 / PI) * cos(cen_lat * 180.0 / PI) / P)) * 180.0 / PI;

Anyhow, It looks like this works.
Code:
double max_lat = -87.5;
double scale   = 1516.17;
double factor  = (double) (width / scale) / (width);
double map_scl = 20.0;
double offset  = 3791.5;
double P       = 0.0;
double C       = 0.0;

Code:
for (int j = 0; j < height; j++)
{

if (j < height / 2)
{
max_lat = max_lat - factor;
}
else
{
max_lat = max_lat + factor;
}

for (int i = 0; i < width; i++)
{
val = (short) (a[i] * 0.5) + rad;
x = (double) (i - offset) * map_scl;
y = (double) (offset - j) * map_scl;
P = sqrt ((x * x) + (y * y));
C = 2 * atan (P / 2 * rad) * 180.0 / PI;
lon = atan2 (x,y) * 180.0 / PI;
// ARCSIN[COS(C) * SIN(LatP) + y * SIN(C) * COS(LatP) / P]
//                lat = asin (cos(C * 180.0 / PI) * sin(cen_lat * 180.0 / PI) + (y * sin(C * 180.0 / PI) * cos(cen_lat * 180.0 / PI) / P)) * 180.0 / PI;

if (i < width / 2)
{
lat = max_lat - (factor * (i - offset));
}
else
{
lat = max_lat + (factor * (i - offset));
}

if (j == 3792)
{
file2 << lat << " " << lon << " " << val << "\n";
}

Not sure why the formula didn't work. I've been all day. cartrite

Top Post subject: Re: LRO - WAC DEM Posted: Wed, 09-10-13, 3:55 GMT  Joined: Fri, 03-04-09, 8:21 GMT
Posts: 221
cartrite wrote:
The top line was the formula from http://pds-geosciences.wustl.edu/lunar/ ... _polar.cat
The bottom line was what I tried. It produced a list of latitudes way out of range.

I tried those formulae out in a spreadsheet first, and found the same issue.
FWIW, here's what I used:
Code:
float x = (sample - projection.SAMPLE_PROJECTION_OFFSET) * projection.MAP_SCALE;
float y = (projection.LINE_PROJECTION_OFFSET - line) * projection.MAP_SCALE;

float p = sqrt( pow(x,2) + pow(y,2) );

latitude = projection.CENTER_LATITUDE + degrees(atan( p / polar_radius ) );

longitude = degrees(atan( x / y ));
if (y < 0)
{
longitude = longitude + 180;
}

Top Post subject: Re: LRO - WAC DEM Posted: Wed, 09-10-13, 10:42 GMT Site Admin Joined: Fri, 31-08-07, 7:01 GMT
Posts: 4675
Location: Hamburg, Germany
As a word of caution:

when inverse trig functions are used (arc tangent,...), one has to be VERY careful about their range, since one deals here with multivalued functions (like the logarithm). Often, using the atan2(y,x) function is preferrable over atan(y/x)!

atan(y/x) returns the principal value of the arc tangent of y/x, expressed in radians, in the interval [-π/2,+π/2].

atan2(y,x) returns the principal value of the arc tangent of y/x, expressed in radians, in the interval [-π,+π].

If the signs of x and y are known individually, one computes y/x via atan2(y,x). It is VERY easy to get required additions of π wrong by using the principal value function atan() from math libraries e.g. in C++. So the true value of an angle is

true arc tan = atan() + n * π

n = ...-1, 0 , 1 , ...

and the user must find out what n is!!

Fridger

Top Post subject: Re: LRO - WAC DEM Posted: Wed, 09-10-13, 11:20 GMT  Joined: Mon, 03-09-07, 23:01 GMT
Posts: 432
Location: Tuscany, Tyrrhenian Sea
There is somethings that I do not understand:

Quote:
C = 2 * atan (P / 2 * rad) * 180.0 / PI;

above C is in degrees; while below the cos and sin does process degrees for then to get degrees once more.

Quote:
lat = asin (cos(C * 180.0 / PI) * sin(cen_lat * 180.0 / PI) + (y * sin(C * 180.0 / PI) * cos(cen_lat * 180.0 / PI) / P)) * 180.0 / PI;

For me it should be:

Code:
C = 2 * atan (P / 2 * rad)

and then:

Code:
lat = asin (cos(C) * sin(cen_lat) + (y * sin(C) * cos(cen_lat) / P)) * 180.0 / PI Top Post subject: Re: LRO - WAC DEM Posted: Wed, 09-10-13, 12:06 GMT  Joined: Thu, 25-10-07, 15:20 GMT
Posts: 992
Location: NE PA, USA
fenerit wrote:
There is somethings that I do not understand:

Quote:
C = 2 * atan (P / 2 * rad) * 180.0 / PI;

above C is in degrees; while below the cos and sin does process degrees for then to get degrees once more.

Quote:
lat = asin (cos(C * 180.0 / PI) * sin(cen_lat * 180.0 / PI) + (y * sin(C * 180.0 / PI) * cos(cen_lat * 180.0 / PI) / P)) * 180.0 / PI;

For me it should be:

Code:
C = 2 * atan (P / 2 * rad)

and then:

Code:
lat = asin (cos(C) * sin(cen_lat) + (y * sin(C) * cos(cen_lat) / P)) * 180.0 / PI Yes, that is what I think it should be too. But it doesn't work so I was trying to see where I was making a mistake (typo) with a missing "(" or a missing ")" or forgetting some conversion.
Quote:
lat = asin (cos(C * 180.0 / PI) * sin(cen_lat * 180.0 / PI) + (y * sin(C * 180.0 / PI) * cos(cen_lat * 180.0 / PI) / P)) * 180.0 / PI;
Doesn't work.
Code:
C = 2 * atan (P / 2 * rad)

and then:

Code:
lat = asin (cos(C) * sin(cen_lat) + (y * sin(C) * cos(cen_lat) / P)) * 180.0 / PI

Doesn't work either.
What does work is calculating the size of a pixel in degrees and adding or subtracting from the maximum latitude depending what pixel it is on the x,y grid.
Code:
..............
double max_lat = -87.5;
double factor  = (double) (width / scale) / (width);
double scale   = 1516.17; // pixels/degree

......................

for (int j = 0; j < height; j++)
{
/*
if (j < height / 2)
{
max_lat = max_lat - factor;
}
else
{
max_lat = max_lat + factor;
}
*/

...............................

for (int i = 0; i < width; i++)

/*
if (i < width / 2)
{
lat = max_lat - (factor * (i - offset));
}
else
{
lat = max_lat + (factor * (i - offset));
}
*/

Still trying to figure where I'm going wrong with
Code:
C = atan (P / (2 * rad)) * 2;
lat = asin(cos((C) * sin(cen_lat)) + (y * sin(C) * cos(cen_lat) / P)) * 180.0 / PI; cartrite

Top Post subject: Re: LRO - WAC DEM Posted: Wed, 09-10-13, 12:51 GMT  Joined: Thu, 25-10-07, 15:20 GMT
Posts: 992
Location: NE PA, USA
After much weeping and gnashing of teeth, also much (like that smiley), The formula for latitude I was trying to write in c++'
Code:
ARCSIN[COS(C) * SIN(LatP) + y * SIN(C) * COS(LatP) / P]

is this:
Code:
P = sqrt ((x * x) + (y * y));
lon = atan2 (x,y) * 180.0 / PI;
C = 2 * atan (P / (2 * rad));
// ARCSIN[COS(C) * SIN(LatP) + y * SIN(C) * COS(LatP) / P]
lat = asin(cos((C) * sin(cen_lat * PI / 180)) + (y * sin(C) * cos(cen_lat * PI / 180) / P)) * 180.0 / PI;

What I was doing wrong was trying to take the sine and cosine from -90 degrees in degrees when I should have used -90 degrees in radians. The value for (cen_lat) is -90.
cartrite

Top Display posts from previous: All posts1 day7 days2 weeks1 month3 months6 months1 year Sort by AuthorPost timeSubject AscendingDescending  Page 2 of 3 [ 33 posts ] Go to page Previous  1, 2, 3  Next

 All times are UTC

#### Who is online

Users browsing this forum: No registered users and 1 guest

 You cannot post new topics in this forumYou cannot reply to topics in this forumYou cannot edit your posts in this forumYou cannot delete your posts in this forumYou cannot post attachments in this forum

Search for:
 Jump to:  Select a forum ------------------ General    The Celestial Matters Lounge    Docs & Tutorials    Artwork, Modelling & Design Development    celestia.Sci Development    Celestial Space    Texture Tools    Cosmological Visualization Project Workshops    Extrasolar Planets    Spice Kernel Files 