Problems recovering 3d mesh from GeoTIFFs created from PDS IMG/LBL files using gdal_translate

Question

Update:

I've discovered that the problem lies in the pixel scale. For example, when I use

`gdalinfo LDEM_512_90S_45S_000_090.tif`

it shows the following for the pixel scale values:

Pixel Size = (59.225293800000003,-59.225293800000003)

However when I use the following code to extract it in my software:

int16_t _count;
double* pixelScale;
TIFFGetField(tif, GTIFF_PIXELSCALE, &_count, &pixelScale);
std::cout << Pixel Scale = (" << pixelScale[0] << ", " << pixelScale[1] <<")\n"

What is printed is:

Pixel Scale = (59.225293800000003, 59.225293800000003)

This is noticeably the right number, but the incorrect sign. When I manually set that value to be negative, the mesh is generated correctly (which makes sense, as a sign flip explains the weird mirroring I was observing).

I'm now investigating why my way of loading the pixel scale (and possibly other values?) is wrong.

---

Original Post

I'm working on converting some lunar data into 3d meshes. I'm using the following set of data: http://imbrium.mit.edu/BROWSE/LOLA_GDR/CYLINDRICAL/ELEVATION/

For each of these data products I download the IMG/LBL pair and convert them to a GeoTIFF via the command line using the following gdal_translate command:

gdal_translate -of GTiff LDEM_4.LBL LDEM_4.tif

I've done this for the LDEM_4, LDEM_128, and all sets of the LDEM_512 DEMs.

The LDEM_4 and LDEM_128 DEMs both give the following modelTiePoint array:

[0, 0, 0, -5.4582e+06, 2.7291e+06, 0]

While here are the mode tie points I've extracted from the LDEM_512 DEMs:

LDEM_512_00N_45N_000_090: [0, 0, 0, -5.4582e+06,  1.36455e+06, 0]
LDEM_512_00N_45N_090_180: [0, 0, 0, -2.7291e+06,  1.36455e+06, 0]
LDEM_512_00N_45N_180_270: [0, 0, 0, -0,           1.36455e+06, 0]
LDEM_512_00N_45N_270_360: [0, 0, 0,  2.7291e+06,  1.36455e+06, 0]
LDEM_512_45N_90N_000_090: [0, 0, 0, -5.4582e+06,  2.7291e+06,  0]
LDEM_512_45N_90N_090_180: [0, 0, 0, -2.7291e+06,  2.7291e+06,  0]
LDEM_512_45N_90N_180_270: [0, 0, 0, -0,           2.7291e+06,  0]
LDEM_512_45N_90N_270_360: [0, 0, 0,  2.7291e+06,  2.7291e+06,  0]
LDEM_512_45S_00S_000_090: [0, 0, 0, -5.4582e+06,  0,           0]
LDEM_512_45S_00S_090_180: [0, 0, 0, -2.7291e+06,  0,           0]
LDEM_512_45S_00S_180_270: [0, 0, 0, -0,           0,           0]
LDEM_512_45S_00S_270_360: [0, 0, 0,  2.7291e+06,  0,           0]
LDEM_512_90S_45S_000_090: [0, 0, 0, -5.4582e+06, -1.36455e+06, 0]
LDEM_512_90S_45S_090_180: [0, 0, 0, -2.7291e+06, -1.36455e+06, 0]
LDEM_512_90S_45S_180_270: [0, 0, 0, -0,          -1.36455e+06, 0]
LDEM_512_90S_45S_270_360: [0, 0, 0,  2.7291e+06, -1.36455e+06, 0]

In my C++ program I'm using libgeotiff and proj to generate a 3d mesh. My general flow is to loop over all pixels (assuming each tif is in row-major order), convert those to model space using the PIXELSCALE and TIEPOINTS, then use proj_trans() to calculate the lat/long values. From there, I use the semi-major and semi-minor axes to calculate the real-world cartesian coordinate. Here is what that code looks like:

// Read the TIF:
TIFF* tif = XTIFFOpen(filepath.c_str(), "r");
GTIF* gtif = GTIFNew(tif);
GTIFDefn* psDefn = GTIFAllocDefn();

// Code to read the height data is long and so ommitted from this post:
std::vector<float> heights = readHeightData(tif);

// Get the model tie point:
int16_t _count;
double* modelTiePoint;
TIFFGetField(tif, GTIFF_TIEPOINTS, &_count, &modelTiePoint);

// Get the pixel scale:
double* pixelScale;
TIFFGetField(tif, GTIFF_PIXELSCALE, &_count, &pixelScale);

// Get the geotiff definition:
GTIFGetDefn(gtif, psDefn);

// Get the source CRS proj4 definition:
char* _pszProjection = GTIFGetProj4Defn(psDefn);
std::string pszProjection = std::string(_pszProjection);
CPLFree(_pszProjection);

// Get the semi-major axis:
double a = psDefn->SemiMajor;
double b = psDefn->SemiMinor;

// Create the target CRS proj4 string:
std::string szLongLat = "+proj=longlat ";
char _szLongLat[512];
sprintf_s(_szLongLat, "+a=%.3f +b=%.3f", a, b);
szLongLat = szLongLat + std::string(_szLongLat);

GTIFFreeDefn(psDefn);
GTIFFree(gtif);
XTIFFClose(tif);

// Setup proj CRS-to-CRS
PJ_CONTEXT* ctx = proj_context_create();
PJ* psPJ = proj_create_crs_to_crs(ctx, pszProjection.c_str(), szLongLat.c_str(), NULL);

// Loop over all pixels:
for (uint32_t i = 0; i < width; i++){
    for (uint32_t j = 0; j < height; j++){
        PJ_COORD coord;
        coord.xyzt.x = (static_cast<double>(i) * pixelScale[0]) + modelTiePoint[3];
        coord.xyzt.y = (static_cast<double>(j) * pixelScale[1]) + modelTiePoint[4];
        coord.xyzt.z = 0;
        coord.xyzt.t = 0;

        coord = proj_trans(psPJ, PJ_FWD, coord);

        double lon = coord.xyzt.x;
        double lat = coord.xyzt.y;
        double alt = static_cast<double>(heights[ i + (j*width) ]);

        positions[idx] = ellipsoidToCartesian(lon, lat, alt, a, b);
    }
}

I should note that the pszProjection I extract is always the same for ALL of the GeoTIFFs (I'm not sure if this is relevant, but I figured I'd point it out). This is what it is:

+proj=eqc +lat_0=0.000000000 +lon_0=180.000000000 +lat_ts=0.000000000 +x_0=0.000 +y_0=0.000 +a=1737400.000 +b=1737400.000 +units=m

In any case, the above code appears to work correctly for both the LDEM_4 and LDEM_128 datasets. They both produce sensible models of the moon and match one another. However LDEM_512 is returning a garbled mess. To illustrate this, I exported the 3d meshes my software has produced into OBJ files that I then imported into Blender to better visualize:

It appears what is happening is that is:

• The 45N_90N_00_90,45N_90N_90_180,45N_90N_180_270,45N_90N_270_360 DEMs appear to be 100% correct (at least, they match exactly with the LDEM_4 and LDEM_128 models I'm producing).
• The 0N_45N_00_90, 0N_45N_90_180, 0N_45N_180_270, 0N_45N_270_360 DEMs are being treated as if they start with a latitude of 45N, and have a longitude 180 degrees off what they should. (You can see in the above image some artifacting towards the top of the LDEM_512 model. That is because the models from these DEMs are sitting directly ontop of the correctly placed 45N_* DEMs).
• The 45S_0S_00_90, 45S_0S_90_180, 45S_0S_180_270, 45S_0S_270_360 DEMs are mirrored over the equator (as if spanning from 0N->45N instead of 45S->0S). Additionally, they also have their longitude offset by 180 degrees.
• The 90S_45S_00_90,90S_45S_90_180,90S_45S_180_270,90S_45S_270_360 DEMs are spanning from 45S->0S instead of from 90S->45S. They ALSO have their longitude offset by 180 degrees from what they should.

I'm not sure where the problem here is. I think it might be with the tie points but I can't quite figure out exactly what is going wrong, or how to fix it.

Answer 1

After examining many cases I noticed that only the y-component of my pixelScale value was wrong, and it only ever had a sign error. Looking through some libgeotiff code I found this line:

...
*x = (*x - tiepoints[0]) * pixel_scale[0] + tiepoints[3];
*y = (*y - tiepoints[1]) * (-1 * pixel_scale[1]) + tiepoints[4];
...

Which implies that the y-component of pixel_scale is to be negated. Making this change has fixed all of my issues. I'm unsure why the value obtained using TIFFGetField is negative, but this resolved my issue completely.