nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Getting the Elvation from a surface raster

JosephSheffield

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

Encuentra más publicaciones etiquetadas con

Accepted answers

All comments

BKuiper

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.

NeilClemmons_GISinc

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

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

Untitled.png

GIS und 3D bei der RAG.PDF

BKuiper

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!

JosephSheffield

Thanks for your help!

NeilClemmons_GISinc

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

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_GISinc

Post the exact code you're using now.

JosephSheffield

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.

NeilClemmons_GISinc

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

JosephSheffield

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

NeilClemmons_GISinc

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.

JosephSheffield

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.

NeilClemmons_GISinc

The example link is at the bottom of the page:

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

JosephSheffield

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.

NeilClemmons_GISinc

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

JosephSheffield

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.

Únete a la Conversación

Los miembros registrados pueden publicar, seguir actualizaciones y más. ¿Nuevo aquí? Regístrate gratis.

Enlaces Rápidos

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.