1

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:

Comparison

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.

1 Answer 1

0

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.

1
  • It's because in pixel space the Y coordinates (row numbers) increase from top -> bottom i.e row 0 is top of the image, but map Y coordinates (Latitude or Northing) increase bottom -> top.
    – user2856
    Jun 4, 2023 at 13:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.