It is currently Tue, 25-07-17, 2:44 GMT

All times are UTC




Post new topic Reply to topic  [ 33 posts ]  Go to page Previous  1, 2, 3  Next
Author Message
 Post subject: Re: LRO - WAC DEM
PostPosted: Sun, 06-10-13, 16:28 GMT 
Offline
User avatar

Joined: Fri, 03-04-09, 8:21 GMT
Posts: 211
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
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Sun, 06-10-13, 17:04 GMT 
Offline
User avatar

Joined: Thu, 25-10-07, 15:20 GMT
Posts: 991
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
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Sun, 06-10-13, 19:21 GMT 
Offline
User avatar

Joined: Tue, 04-09-07, 21:55 GMT
Posts: 757
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
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Tue, 08-10-13, 13:01 GMT 
Offline
User avatar

Joined: Fri, 03-04-09, 8:21 GMT
Posts: 211
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
000 - south pole 01.jpg [ 391.17 KiB | Viewed 2049 times ]

Attachment:
001 - terrain 01.jpg
001 - terrain 01.jpg [ 199.62 KiB | Viewed 2049 times ]

Attachment:
001 - terrain 03.jpg
001 - terrain 03.jpg [ 107.79 KiB | Viewed 2049 times ]


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

Top
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Tue, 08-10-13, 13:07 GMT 
Offline
User avatar

Joined: Fri, 03-04-09, 8:21 GMT
Posts: 211
And some additional views (CLICK for fullsize)....

Shackleton Base up close:
Attachment:
002 - terrain 03.jpg
002 - terrain 03.jpg [ 136.03 KiB | Viewed 2049 times ]


A view from high up inside the dome overlooking buildings far below and the lunar terrain outside:
Attachment:
004 - Inside Looking Out.jpg
004 - Inside Looking Out.jpg [ 121.35 KiB | Viewed 2047 times ]


Shackleton Base and terrain as it appears from the Malapert Mountains communication centre:
Attachment:
003 - from malapert.jpg
003 - from malapert.jpg [ 88.32 KiB | Viewed 2049 times ]


Top
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Tue, 08-10-13, 16:28 GMT 
Offline
User avatar

Joined: Thu, 25-10-07, 15:20 GMT
Posts: 991
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
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Wed, 09-10-13, 1:03 GMT 
Offline
User avatar

Joined: Fri, 03-04-09, 8:21 GMT
Posts: 211
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
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Wed, 09-10-13, 1:28 GMT 
Offline
User avatar

Joined: Thu, 25-10-07, 15:20 GMT
Posts: 991
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
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Wed, 09-10-13, 1:41 GMT 
Offline
Site Admin
User avatar

Joined: Fri, 31-08-07, 7:01 GMT
Posts: 4467
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

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
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Wed, 09-10-13, 2:25 GMT 
Offline
User avatar

Joined: Thu, 25-10-07, 15:20 GMT
Posts: 991
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. :shock:
cartrite


Top
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Wed, 09-10-13, 3:55 GMT 
Offline
User avatar

Joined: Fri, 03-04-09, 8:21 GMT
Posts: 211
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
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Wed, 09-10-13, 10:42 GMT 
Offline
Site Admin
User avatar

Joined: Fri, 31-08-07, 7:01 GMT
Posts: 4467
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
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Wed, 09-10-13, 11:20 GMT 
Offline
User avatar

Joined: Mon, 03-09-07, 23:01 GMT
Posts: 385
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


:rool:


Top
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Wed, 09-10-13, 12:06 GMT 
Offline
User avatar

Joined: Thu, 25-10-07, 15:20 GMT
Posts: 991
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


:rool:

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;

:wall:
cartrite


Top
 Profile  
 
 Post subject: Re: LRO - WAC DEM
PostPosted: Wed, 09-10-13, 12:51 GMT 
Offline
User avatar

Joined: Thu, 25-10-07, 15:20 GMT
Posts: 991
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
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 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 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

Search for:
Jump to:  
cron
Powered by phpBB® Forum Software © phpBB Group