• View topic - LRO

View unanswered posts | View active topics

Board index » Development » Celestial Space

All times are UTC

LRO - WAC DEM

Moderator: Moderators

Page 2 of 3

[ 33 posts ]

Go to page Previous 1, 2, 3 Next

Print view

Previous topic | Next topic

Author

Message

chuft-captain

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

cartrite

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

John Van Vliet

Post subject: Re: LRO - WAC DEM

Posted: Sun, 06-10-13, 19:21 GMT

Joined: Tue, 04-09-07, 21:55 GMT
Posts: 859
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

ad a burn mask to the image , add the two images together
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

chuft-captain

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:

000 - south pole 01.jpg [ 391.17 KiB | Viewed 5633 times ]

Attachment:

001 - terrain 01.jpg [ 199.62 KiB | Viewed 5633 times ]

Attachment:

001 - terrain 03.jpg [ 107.79 KiB | Viewed 5633 times ]

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

Top

chuft-captain

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:

002 - terrain 03.jpg [ 136.03 KiB | Viewed 5633 times ]

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

Attachment:

004 - Inside Looking Out.jpg [ 121.35 KiB | Viewed 5631 times ]

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

Attachment:

003 - from malapert.jpg [ 88.32 KiB | Viewed 5633 times ]

Top

cartrite

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

chuft-captain

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

cartrite

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

t00fri

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: 4666
Location: Hamburg, Germany

cartrite wrote:

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

angle_radian = angle_deg * Pi/180.0
angle_deg = angle_radian * 180.0/Pi

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

Good luck,
Fridger

Top

cartrite

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 rad     = 1737400.0;
    double offset  = 3791.5;
    double P       = 0.0;
    double C       = 0.0;

Code:

        for (int j = 0; j < height; j++)
       {
            fread(a, 2, width, file1);

            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 :wall:

all day.

cartrite

Top

chuft-captain

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

t00fri

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: 4666
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

fenerit

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

cartrite

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++)
       {
            fread(a, 2, width, file1);
/*
            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

cartrite

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 :wall:

(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

Page 2 of 3

[ 33 posts ]

Go to page Previous 1, 2, 3 Next

Board index » Development » Celestial Space

All times are UTC

Who is online

Users browsing this forum: No registered users and 2 guests

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