Getting the Elvation from a surface raster

JosephSheffield · ‎10-22-2010

I'm trying to get the elevation value from a raster along a specific line. I'm passing in the raster layer as a IRasterLayer and the feature class as IFeatureLayer. I know I have to the raster from a IRasterLayer data type to an ISurface data type. Does anyone know the way to do so?

Here is the relevant code block:

public void GetThicknessData(IRasterLayer rAlamoSurface, IFeatureLayer fThicknessProfile)
          {
              try
              {
                  IFeature fProfile = fThicknessProfile.FeatureClass.GetFeature(11);
                  IPointCollection pPointCollection = null;
                  IPolyline pPolyLine = null;
                  IGeometry pGeom = null;
                  Double dCellSize = 30;
                  IRasterSurface pRasterSurface = rAlamoSurface as IRasterSurface;

                  ISurface pSurface = pRasterSurface as ISurface;
                  pPolyLine = fProfile.Shape as IPolyline;
                  pPointCollection = pPolyLine as IPointCollection;
                  pGeom = pPointCollection as IGeometry;

                  IPoint tmpPt = null;
                  IPoint tmpPt2 = null;
                  IPoint stepPoint = null;
                  int ArrayCtr, loopCount, ptCtr, sgCtr;

                  ArrayCtr = 0;
                  double[] XArray = new double[1000];
                  double[] YArray = new double[1000];
                  double[] ElevArray = new double[1000];
                  double pElev, dX, dY, segLength, runStep, riseStep, lastPointX, lastPointY, gridNoData;
                  gridNoData = -99999;

                  tmpPt = pPointCollection.get_Point(0);
                  XArray[ArrayCtr] = tmpPt.X;
                  YArray[ArrayCtr] = tmpPt.Y;

                  pElev = pSurface.GetElevation(tmpPt);
                  if (pElev != gridNoData)
                  {
                      ElevArray[ArrayCtr] = pElev;
                  }

                  for (ptCtr = 1; ptCtr <= pPointCollection.PointCount - 1; ptCtr++)
                  {
                      tmpPt.X = pPointCollection.get_Point(ptCtr - 1).X;
                      tmpPt.Y = pPointCollection.get_Point(ptCtr - 1).Y;
                      tmpPt2.X = pPointCollection.get_Point(ptCtr).X;
                      tmpPt2.Y = pPointCollection.get_Point(ptCtr).Y;

                      dX = (tmpPt2.X - tmpPt.X);
                      dY = (tmpPt2.Y - tmpPt.Y);

                      segLength = Math.Sqrt(Math.Pow(dX, 2) + Math.Pow(dY, 2));
                      runStep = dCellSize * (dX) / segLength;
                      riseStep = dCellSize * (dY) / segLength;

                      loopCount = (int)Math.Round((segLength * 10000) / (dCellSize * 10000));
                      //*10000 because the cell size is so small it causes a divide by zero error.
                      lastPointX = tmpPt.X;
                      lastPointY = tmpPt.Y;

                      for (sgCtr = 1; sgCtr <= loopCount; sgCtr++)
                      {
                          ArrayCtr = ArrayCtr + 1;
                          stepPoint.X = lastPointX + runStep;
                          stepPoint.Y = lastPointY + riseStep;
                          XArray[ArrayCtr] = stepPoint.X;
                          YArray[ArrayCtr] = stepPoint.Y;
                          pElev = pSurface.GetElevation(stepPoint);

                          if (pElev != gridNoData)
                          {
                              ElevArray[ArrayCtr] = pElev;
                              lastPointX = stepPoint.X;
                              lastPointY = stepPoint.Y;
                          }
                      }

                      if ((segLength * 10000) % (dCellSize * 10000) != 0)
                      {
                          ArrayCtr = ArrayCtr + 1;
                          stepPoint.X = lastPointX + runStep;
                          stepPoint.Y = lastPointY + riseStep;
                          XArray[ArrayCtr] = stepPoint.X;
                          YArray[ArrayCtr] = stepPoint.Y;
                          pElev = pSurface.GetElevation(stepPoint);
                          if (pElev != gridNoData)
                          {
                              ElevArray[ArrayCtr] = pElev;
                          }
                      }
                  }
              }

              catch (Exception ex)
              {
                  MessageBox.Show(ex.Message);
              }
             }

NeilClemmons · ‎11-05-2010

Post the exact code you're using now.

JosephSheffield · ‎11-05-2010

public void GetThicknessData(IRasterLayer rAlamoSurface, IFeatureLayer fThicknessProfile)
{
try
{
IFeature fProfile = fThicknessProfile.FeatureClass.GetFeature(11);
IPointCollection pPointCollection = null;
IPolyline pPolyLine = null;
IGeometry pGeom = null;
Double dCellSize = 30;
IRasterSurface pRasterSurface = null;
pRasterSurface.RasterBand = rasterBand;

ISurface pSurface = pRasterSurface as ISurface;
pPolyLine = fProfile.Shape as IPolyline;
pPointCollection = pPolyLine as IPointCollection;
pGeom = pPointCollection as IGeometry;

IPoint tmpPt = null;
IPoint tmpPt2 = null;
IPoint stepPoint = null;
int ArrayCtr, loopCount, ptCtr, sgCtr;

ArrayCtr = 0;
double[] XArray = new double[1000];
double[] YArray = new double[1000];
double[] ElevArray = new double[1000];
double pElev, dX, dY, segLength, runStep, riseStep, lastPointX, lastPointY, gridNoData;
gridNoData = -99999;

tmpPt = pPointCollection.get_Point(0);
XArray[ArrayCtr] = tmpPt.X;
YArray[ArrayCtr] = tmpPt.Y;

pElev = pSurface.GetElevation(tmpPt);
if (pElev != gridNoData)
{
ElevArray[ArrayCtr] = pElev;
}

for (ptCtr = 1; ptCtr <= pPointCollection.PointCount - 1; ptCtr++)
{
tmpPt.X = pPointCollection.get_Point(ptCtr - 1).X;
tmpPt.Y = pPointCollection.get_Point(ptCtr - 1).Y;
tmpPt2.X = pPointCollection.get_Point(ptCtr).X;
tmpPt2.Y = pPointCollection.get_Point(ptCtr).Y;

dX = (tmpPt2.X - tmpPt.X);
dY = (tmpPt2.Y - tmpPt.Y);


segLength = Math.Sqrt(Math.Pow(dX, 2) + Math.Pow(dY, 2));
runStep = dCellSize * (dX) / segLength;
riseStep = dCellSize * (dY) / segLength;

loopCount = (int)Math.Round((segLength * 10000) / (dCellSize * 10000));
//*10000 because the cell size is so small it causes a divide by zero error.
lastPointX = tmpPt.X;
lastPointY = tmpPt.Y;

for (sgCtr = 1; sgCtr <= loopCount; sgCtr++)
{
ArrayCtr = ArrayCtr + 1;
stepPoint.X = lastPointX + runStep;
stepPoint.Y = lastPointY + riseStep;
XArray[ArrayCtr] = stepPoint.X;
YArray[ArrayCtr] = stepPoint.Y;
pElev = pSurface.GetElevation(stepPoint);

if (pElev != gridNoData)
{
ElevArray[ArrayCtr] = pElev;
lastPointX = stepPoint.X;
lastPointY = stepPoint.Y;
}
}

if ((segLength * 10000) % (dCellSize * 10000) != 0)
{
ArrayCtr = ArrayCtr + 1;
stepPoint.X = lastPointX + runStep;
stepPoint.Y = lastPointY + riseStep;
XArray[ArrayCtr] = stepPoint.X;
YArray[ArrayCtr] = stepPoint.Y;
pElev = pSurface.GetElevation(stepPoint);
if (pElev != gridNoData)
{
ElevArray[ArrayCtr] = pElev;
}
}
}
}


catch (Exception ex)
{
MessageBox.Show(ex.Message);
}
}

NeilClemmons · ‎11-05-2010

IRasterSurface pRasterSurface = null;
pRasterSurface.RasterBand = rasterBand;

ISurface pSurface = pRasterSurface as ISurface;

You're trying to assign a value (rasterBand) to a property on an object (pRasterSurface) that hasn't been instantiated. It's value is null so that's why you're getting that error. Also, where is rasterBand declared and where is it being set? What you should be doing is assigning the value returned by GetSurfaceFromLayer to pSurface:

//IRasterSurface pRasterSurface = null; delete this
//pRasterSurface.RasterBand = rasterBand; delete this
ISurface pSurface = GetSurfaceFromLayer(rAlamoSurface);

Get rid of pRasterSurface altogether. You'll also need to translate the GetSurfaceFromLayer method into C# and include it in your code project.

JosephSheffield · ‎11-05-2010

Thanks for your help!

BKuiper · ‎03-15-2012

I'm using this code to replace an existing module that uses GlobalMapper software, but the Esri code is much much slower. Are there any other ways to get elevation data from a raster that support interpolation but is much quicker.

Thank you!

BKuiper · ‎03-16-2012

Here is a screenshot of my proof of concept.

The first one shows the Esri function, the second one the GlobalMapper function.

I'm trying to get the elevation for 10.000 data points as quickly as possible.

Is there any way i can improve the speed using the Esri library? this is just a pathetic performance.

So currently i'm using the Esri function ISurface.GetElevation()

NeilClemmons · ‎03-16-2012

It may be possible to speed things up...

If you're using ArcObjects in a .NET application it will be inherently slow because .NET does not natively support COM. In order to use ArcObjects you have to go through COM Interop which will slow down any ArcObjects method call (not to mention that .NET is just slower in general than some other language platforms). One solution is to write the code that needs to execute quickly in a language that natively supports COM (such as C++ or VB6) and compile that into a dll that you can call from your .NET application. Another solution is to code the needed algorithms yourself. ESRI supports a wide variety of datasources, platforms, etc. so their code cannot be written for a specific environment and therefore may not run as efficiently as it could if it were written for something specific. Depending on what datasource you're using, you may be able to write your own code that takes advantage of shortcuts specific to that data format (or find code on the internet).

BKuiper · ‎03-16-2012

Hi Neil,

thanks for taking the time to respond on my question.

I understand what you are saying. And writing my own library is always an option.

But in this case i'm also using a COM wrapped call to Global Mapper in .NET which is about 2000 (!) times faster then Esri ArcObjects. I would at least like to get close to the same performance as GlobalMapper using Esri calls.

I'm hoping this is possible.

Currently i'm using IFunctionalSurface2.get_Z(lon, lat) to get the Z, but i'm looking to each DEM separately. I think I could probably get some performance gain if i could merge the rasters (in memory, during runtime) and only have to call .get_Z once on the combined raster (ISurface) instead of calling it on each raster and see if it was within the extent of the raster.