topic Re: Getting the Elvation from a surface raster in ArcObjects SDK Questions

Getting the Elvation from a surface raster

JosephSheffield — Fri, 22 Oct 2010 13:02:48 GMT

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);
              }
             }

Re: Getting the Elvation from a surface raster

JosephSheffield — Sat, 11 Dec 2021 10:04:58 GMT

So I've figure out part of the problem. I can get the raster to a IRasterBand class, but I can't figure out the code/syntax to get the raster to a IRasterSurface class: Here is the code:

IFeature fProfile = fThicknessProfile.FeatureClass.GetFeature(11);
                  IPointCollection pPointCollection = null;
                  IPolyline pPolyLine = null;
                  IGeometry pGeom = null;
                  Double dCellSize = 30;
                  IRasterBandCollection pRasterBands = null;
                  pRasterBands = rAlamoSurface.Raster as IRasterBandCollection;
                  IRasterBand pRasterBand = pRasterBands.Item(0);
                  IRasterSurface pRasterSurface = null;
                  //I get the error here....trying to cast an IRasterBand to an IRasterSurface
                  pRasterSurface.RasterBand = pRasterBand as IRasterSurface;
                  //I've found other VB.net references that use:
                  //pRasterSurface.RasterBand = pRasterBand

Thanks in advance for anyone who can help me.

Re: Getting the Elvation from a surface raster

NeilClemmons — Mon, 25 Oct 2010 13:14:59 GMT

If you open the developer help to the ISurface.GetElevation method topic you'll find a link to the Get Elevation developer sample.

Re: Getting the Elvation from a surface raster

JosephSheffield — Fri, 05 Nov 2010 01:01:25 GMT

The Developer Help does not include sample code that helps me figure out how to cast the IRasterBand to IRasterSurface. I have the code working in VBA for Arc 9.3, I'm migrating it over to Arc 10. Every resource I can finds leads me to believe that the correct syntax is:

pRasterSurface.RasterBand = pRasterBand

However, when I run the code nothing is passed into pRasterSurface. I get a null value.

Re: Getting the Elvation from a surface raster

NeilClemmons — Fri, 05 Nov 2010 12:17:05 GMT

The example link is at the bottom of the page:

http://resources.esri.com/help/9.3/arcgisserver/apis/arcobjects/esrigeodatabase/ISurface_GetElevation.htm

Re: Getting the Elvation from a surface raster

JosephSheffield — Fri, 05 Nov 2010 17:57:50 GMT

I'm using ArcGIS 10. I've used the same methodology as with the code I have for VBA in ArcGIS 9.3.

Based on the link I've tired:

pRasterSurface.RasterBand = pRasterBands.Item(0);

And I still get an error.

Re: Getting the Elvation from a surface raster

NeilClemmons — Fri, 05 Nov 2010 18:02:36 GMT

You don't need to use IRasterSurface at all. All you need is an ISurface reference so that you can call ISurface.GetElevation(). The example in the link above shows how to get the ISurface reference from the IRasterLayer reference.

Re: Getting the Elvation from a surface raster

JosephSheffield — Fri, 05 Nov 2010 19:04:25 GMT

In the same code I see two ways of doing it:

pSurface = GetSurfaceFromLayer(pLayer.Name;

pRasterSurface.RasterBand = pRasterBands.Item(0);

In both cases I get an error. I don't see a method (GetSurfaceFromLayer) on the pSurface.

The error I am getting is "Object reference not set to an instance of an object."

Re: Getting the Elvation from a surface raster

NeilClemmons — Sat, 11 Dec 2021 10:05:01 GMT

GetSurfaceFromLayer is not a method on ISurface; it is a function provided in the example that will return an ISurface reference given a layer name. Here is the same function translated to VB.NET. Pass in a reference to the raster layer and it will return the surface.

        ''' <summary>
        ''' Returns the surface from the layer with the specified name.
        ''' </summary>
        ''' <param name="layer">The layer whose surface will be returned</param>
        ''' <returns></returns>
        ''' <remarks></remarks>
        Public Shared Function GetSurfaceFromLayer(ByVal layer As ILayer) As ISurface
            ' Get the actual surface from the layer.
            If TypeOf layer Is ITinLayer Then
                ' QI
                Dim tinLayer As ITinLayer = CType(layer, ITinLayer)

                ' Return surface from TIN layer.
                Return CType(tinLayer.Dataset, ISurface)
            ElseIf TypeOf layer Is IRasterLayer Then
                ' QI
                Dim rasterLayer As IRasterLayer = CType(layer, IRasterLayer)

                ' Get surface from 3d extension on layer.
                Dim layerExtensions As ILayerExtensions = CType(rasterLayer, ILayerExtensions)
                For i As Integer = 0 To layerExtensions.ExtensionCount - 1
                    Dim dddProperties As I3DProperties = CType(layerExtensions.Extension(i), I3DProperties)
                    If Not dddProperties Is Nothing Then
                        If dddProperties.BaseOption = esriBaseOption.esriBaseSurface Then
                            If Not dddProperties.BaseSurface Is Nothing Then
                                Return CType(dddProperties.BaseSurface, ISurface)
                            End If
                        End If
                        Exit For
                    End If
                Next i

                ' Get surface from raster band collection.
                Dim rasterBandCollection As IRasterBandCollection = CType(rasterLayer.Raster, IRasterBandCollection)
                Dim rasterBand As IRasterBand = rasterBandCollection.Item(0)
                Dim rasterSurface As IRasterSurface = New RasterSurface
                rasterSurface.RasterBand = rasterBand

                Return CType(rasterSurface, ISurface)
            Else
                Return Nothing
            End If
        End Function

Re: Getting the Elvation from a surface raster

JosephSheffield — Fri, 05 Nov 2010 19:42:36 GMT

I see that now.

However, the step that gives me an error each time is:

 rasterSurface.RasterBand = rasterBand;

In every example of code I can find I see this same step and I get the same error each time.

Re: Getting the Elvation from a surface raster

NeilClemmons — Fri, 05 Nov 2010 19:45:27 GMT

Post the exact code you're using now.

Re: Getting the Elvation from a surface raster

JosephSheffield — Sat, 11 Dec 2021 10:05:04 GMT

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);
}
}

Re: Getting the Elvation from a surface raster

NeilClemmons — Fri, 05 Nov 2010 20:56:55 GMT

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.

Re: Getting the Elvation from a surface raster

JosephSheffield — Fri, 05 Nov 2010 22:24:24 GMT

Thanks for your help!

Re: Getting the Elvation from a surface raster

BKuiper — Fri, 16 Mar 2012 03:14:42 GMT

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!

Re: Getting the Elvation from a surface raster

BKuiper — Fri, 16 Mar 2012 12:12:11 GMT

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()

Re: Getting the Elvation from a surface raster

NeilClemmons — Fri, 16 Mar 2012 12:30:06 GMT

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).

Re: Getting the Elvation from a surface raster

BKuiper — Fri, 16 Mar 2012 14:53:00 GMT

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.