http://forum.celestialmatters.org/ 

LRO  WAC DEM http://forum.celestialmatters.org/viewtopic.php?f=9&t=418 
Page 2 of 3 
Author:  chuftcaptain [ Sun, 061013, 16:28 GMT ] 
Post subject:  Re: LRO  WAC DEM 
cartrite wrote: I got this output. value = 6479 Descaled Value = 1734160.5 note* not sure what Descaled Value is? If you notice the line in your gdalinfo output:Possible first line of tab file. The lat/long coords are probably correct. The value still questionable. Code: 86.46504 45.00000 1730.921 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 ) ... in other words, 3239.5 m below the nominal surface.= 1737400 + (6479 * 0.5) = 1737400  3239.5 = 1734160.5 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. 
Author:  cartrite [ Sun, 061013, 17:04 GMT ] 
Post subject:  Re: LRO  WAC DEM 
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://pdsgeosciences.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 
Author:  John Van Vliet [ Sun, 061013, 19:21 GMT ] 
Post subject:  Re: LRO  WAC DEM 
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 
Author:  chuftcaptain [ Tue, 081013, 13:01 GMT ] 
Post subject:  Re: LRO  WAC DEM 
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:

Author:  chuftcaptain [ Tue, 081013, 13:07 GMT ] 
Post subject:  Re: LRO  WAC DEM 
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:

Author:  cartrite [ Tue, 081013, 16:28 GMT ] 
Post subject:  Re: LRO  WAC DEM 
chuftcaptain 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 
Author:  chuftcaptain [ Wed, 091013, 1:03 GMT ] 
Post subject:  Re: LRO  WAC DEM 
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. 
Author:  cartrite [ Wed, 091013, 1:28 GMT ] 
Post subject:  Re: LRO  WAC DEM 
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 
Author:  t00fri [ Wed, 091013, 1:41 GMT ] 
Post subject:  Re: LRO  WAC DEM 
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 
Author:  cartrite [ Wed, 091013, 2:25 GMT ] 
Post subject:  Re: LRO  WAC DEM 
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://pdsgeosciences.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 all day. cartrite 
Author:  chuftcaptain [ Wed, 091013, 3:55 GMT ] 
Post subject:  Re: LRO  WAC DEM 
cartrite wrote: The top line was the formula from http://pdsgeosciences.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; } 
Author:  t00fri [ Wed, 091013, 10:42 GMT ] 
Post subject:  Re: LRO  WAC DEM 
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 
Author:  fenerit [ Wed, 091013, 11:20 GMT ] 
Post subject:  Re: LRO  WAC DEM 
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 
Author:  cartrite [ Wed, 091013, 12:06 GMT ] 
Post subject:  Re: LRO  WAC DEM 
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 
Author:  cartrite [ Wed, 091013, 12:51 GMT ] 
Post subject:  Re: LRO  WAC DEM 
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 
Page 2 of 3  All times are UTC 
Powered by phpBB® Forum Software © phpBB Group http://www.phpbb.com/ 