6

I have to create a terrain heightmap(raster tif preferably) from mesh of xyz ungridded points(in geographical coordinates), I tried taking min and max of bbox of mesh and making width & height by subtracting max-min and then querying for height on all the vertices in the mesh. But it works only if they are in geographical cordinate system n not lat-long as then it can be really tiny. Any ideas how to do it?

My code snippet to create Raster :

double startX = mesh.verticesBBox.bbmin.X;
double startY = mesh.verticesBBox.bbmin.Y;
double bottomX = mesh.verticesBBox.bbmax.X;
double bottomY = mesh.verticesBBox.bbmax.Y;

int width = fabs(bottomX - startX);
int height = fabs(bottomY - startY);

const int channel = 1;  //Dem
const char *pszFormat = "GTiff";
GDALDriver *poDriver;
char **papszMetadata;
poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
GDALDataset *poDstDs;
int count = 0, a =0 , b=0;
poDstDs = poDriver->Create(tfile.toString().c_str(), width, height, channel, GDT_Float32, NULL);
GDALRasterBand* band = poDstDs->GetRasterBand(1);
std::vector<unsigned char> data(width*height*GDALGetDataTypeSize(GDT_Float32) / 8, 0);
int bufferSizeX = width*GDALGetDataTypeSize(GDT_Float32) / 8;
int bufferSizeY = height*GDALGetDataTypeSize(GDT_Float32) / 8;

for (///////////vertex array iterator)
{
    double h;
    dem.GetHeight(Vector2<REAL>(ver->X, ver->Y), h);
    data[count] = unsigned char(h);
    int locationX = fabs(ver->X - startX);
    int locationY = fabs(ver->Y - startY);
    if ((width > locationX) && (height > locationY))
    {
            CPLErr e = band->RasterIO(GF_Write, locationX, locationY, 1, 1, &data[count], 1, 1, GDT_Float32, 0, 0);
    }
        count++;
}
  • XYZ is a supported GDAL format, you should be able to use GDAL_Translate to convert it to a GeoTiff. Provided the points are regularly spaced. If not you will need to convert it to a triangulation of some sort. – Michael Stimson Nov 12 '14 at 22:35
  • Gdal doesnt support ungridded xyz files thatz why have to make a workaround – saagn Nov 12 '14 at 22:36
  • I already have the triangulation but now i am struck what to do next. I tried taking min and max of bbox of mesh and making width & height by subtracting min out of max for x n y but it works only if they(x & y) are in geographical coordinate system n not lat-long as then it can be really tiny. Then i query for height on vertices of the mesh. Any other ideas ? – saagn Nov 12 '14 at 22:49
  • What format is the triangulation in? Could you project the points to a projected coordinate system and re-make the mesh. – Michael Stimson Nov 12 '14 at 22:52
  • Its a custom struct of vector3<doubles>. Yes i can repoject them. – saagn Nov 12 '14 at 22:54
2

Here is some of my C++ code for working with triangles:

Structures, in order to understand my code understanding these structures is necessary:

struct CoPair // linked list of coordinates
{
    long ID;
    float X;
    float Y;
    float Z;
    CoPair *NxtPt;
};

struct Tri // linked list of triangles made from 3 CoPair
{
    bool Active;
    CoPair *A;
    CoPair *B;
    CoPair *C;
    CP Cent;
    Tri *Next;
};

This one is handy and really fast for determining if a given point falls inside a triangle, it doesn't help with the elevation as the point can be anywhere in the triangle.

bool PointInTriangle(CoPair *p, Tri * T)
{
    // using vector product... if AB.AN > 0 && BC.BN > 0 && CA.CN > 0
    // either all true or all false..
    bool RightAB,RightBC,RightCA;
    //DotProd = (DxA * DyC) - (DxC * DyA);

    RightAB = ((T->B->X - T->A->X)*(p->Y - T->A->Y)) - ((p->X - T->A->X)*(T->B->Y - T->A->Y)) <= 0.0f;
    RightBC = ((T->C->X - T->B->X)*(p->Y - T->B->Y)) - ((p->X - T->B->X)*(T->C->Y - T->B->Y)) <= 0.0f;
    RightCA = ((T->A->X - T->C->X)*(p->Y - T->C->Y)) - ((p->X - T->C->X)*(T->A->Y - T->C->Y)) <= 0.0f;
    return (RightAB == RightBC && RightBC == RightCA);
}

To determine the 'z' on the triangle using Barycentric Coordinates:

float findZonTri(CoPair *p, Tri * T)
{
    float Den = ((T->B->Y - T->C->Y)*(T->A->X - T->C->X) + (T->C->X - T->B->X)*(T->A->Y - T->C->Y));
    float LamA = ((T->B->Y - T->C->Y)*(p->X - T->C->X) + (T->C->X - T->B->X)*(p->Y - T->C->Y)) / Den;
    float LamB = ((T->C->Y - T->A->Y)*(p->X - T->C->X) + (T->A->X - T->C->X)*(p->Y - T->C->Y)) / Den;
    float LamC = 1.0f - LamA - LamB;
    float Z = T->A->Z * LamA + T->B->Z * LamB + T->C->Z * LamC;
    return Z;
}

It is important that the point is determined to be inside the triangle before calling findZonTri or the answer is likely to be garbage - also there's much more calculation involved so will slow the process down.

For each point (cell centre) in the raster I iterate through all the triangles and determine if the point is inside then exit that loop and find the Z until all the cells are checked. If all triangles return false then you're not in any triangles so use a NoData value.

Creating a raster using GDAL:

GDALAllRegister();
unsigned long Rows = (Ymax - Ymin) / CellSize;
unsigned long Cols = (Xmax - Xmin) / CellSize;

GDALDriver* IMGdriver = GetGDALDriverManager()->GetDriverByName(DRVname); // get the ERDAS Imagine driver
double GeoTransform[6];
GeoTransform[0] = Xmin;     // Upper Left X
GeoTransform[1] = CellSize; // W-E pixel size
GeoTransform[2] = 0;        // Rotation, 0 if 'North Up'
GeoTransform[3] = Ymax;     // Upper Left Y
GeoTransform[4] = 0;        // Rotation, 0 if 'North Up'
GeoTransform[5] = -CellSize;    // N-S pixel size

// create the dataset
AccumDS = (GDALDataset*) IMGdriver->Create(AccumName,Cols,Rows,1,GDALDataType::GDT_UInt32,NULL);
// put it in the right place in the world
AccumDS->SetGeoTransform(GeoTransform);
AccumDS->SetProjection(SR_Text); // optional, apply the spatial reference so it's not just numbers

I am using a constant cell size, in this case the x-cell and y-cell size is the same.

I hope that helps.

  • @Micheal I had almost similar structure. But how do you calculate the pixeldensity ? My code snippets: 'code' double startX = mesh.verticesBBox.bbmin.X; double startY = mesh.verticesBBox.bbmin.Y; double bottomX = mesh.verticesBBox.bbmax.X; double bottomY = mesh.verticesBBox.bbmax.Y; int width = fabs(bottomX - startX); int height = fabs(bottomY - startY); GDALDataset poDstDs; int count = 0, a =0 , b=0; poDstDs = poDriver->Create(tfile.toString().c_str(), width, height, channel, GDT_Float32, NULL); GDALRasterBand band = poDstDs->GetRasterBand(1); – saagn Nov 13 '14 at 3:40
  • 1
    You need to decide if you want a regular cell size or a specific number of pixels. If you want 100 x 100 cells then width = 100 and height = 100, xCellSize = (MaxX - MinX) / 100 and yCellSize = (MaxY - MinY) / 100. If you want a particular cell size then width = (int)((MaxX - MinX) / xCellSize) and height = (int)((MaxY - MinY) / yCellSize). For creating the dataset width and height are the number of cells in the raster not the geographic extents. – Michael Stimson Nov 13 '14 at 3:47
  • @Micheal I added snippet in question description, i believe alot more readable :) – saagn Nov 13 '14 at 4:10
  • If you're in lat/lon your width will be near to 0. You can't create a dataset with 0 pixels. You need to convert you width and height to cells. See in my code where I convert width/height into rows and cols. The actual width/height of the cells is governed by the GeoTransform which sets the width/height of the cell. You will also find it quicker to write a block or row of cells rather than one at a time. – Michael Stimson Nov 13 '14 at 4:14
  • For example my current cordinates are: when in geographical Corner Coordinates: Upper Left ( -652500.100, -980005.000) Lower Left ( -652500.100,-1045541.000) Upper Right ( -586964.100, -980005.000) Lower Right ( -586964.100,-1045541.000) when in lat long min 15.567351301100075, 50.735922395701373 max 15.605854709902815, 50.756433378887316 so what will be the be the with in both scenarios? n what will be the cellsize ? if it has like 40000 points(vertices) which make 1198095 triangles. – saagn Nov 13 '14 at 16:31
1

GDAL comes with a gdal_grid utility to create a regular grid from scattered data. I'd start with this one since it seems that you are already using GDAL tools.

If you have dense LiDAR data, check out LAStools, such as las2dem or lasgrid.

  • I tried with gdal_grid utility but seems like the actuall gemometry of DEM changes. I tried almost all the interpolation methods. Might be i am doing something wrong there. – saagn Nov 13 '14 at 0:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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