# Getting raster cell indexes from point coordinates

3183
3
10-29-2013 10:55 AM
Highlighted
New Contributor III
Hi,

I have a point shape file (IFeatureLayer) and a raster file (IRaster) taken from an ArcMap project. My goal is to convert from the point coordinates to raster cell indexes. So, the following steps may be necessary

1) get coordinates (XF,YF) of each point in the feature layer
2) convert coordinates (XF,YF) to raster spatial reference (XR,YR)
3) convert from coordinates (XR,YR) to grid cell indexes (kXR,kYR) in raster file

I am not sure how to code these steps by using the ArcObjets .NET. I am using ArcMap 10.1 and coding in C#. If you could send me some hints or snippets for these steps, it would be greatly appreciated.

Tags (5)
3 Replies
Highlighted
Occasional Contributor
pRaster as IRaster2

x = pFeat.Extent.XMax() 'map coordinate
y = pFeat.Extent.YMax()

i = pRaster.ToPixelColumn(x) 'raster coordinate
j = pRaster.ToPixelRow(y)

value = pRaster.GetPixelValue(0, i, j)
Highlighted
New Contributor III
Hi again,

After the help above from Martin, the help from Jamie P. at ESRI Support Team, and a trial and error process, please find below a working C# code that solves this problem. If you find a better way to accomplish this, please let me know. Thanks!

private int getPoints(IRaster raster, IFeatureLayer featureLayer,
out long[] kxR, out long[] kyR)
{
IRasterProps rasterProps = (IRasterProps)raster; // Raster properties

// get coordinates (xF,yF) of each point in the feature layer ---------------------------------------
IQueryFilter queryFilter = new QueryFilterClass(); // Set up the query by default
int nP = featureLayer.FeatureClass.FeatureCount(queryFilter); // number of points
Debug.Print("Number of points in shape file= " + nP.ToString());

double[] xF = new double[nP];
double[] yF = new double[nP];
for (int kP = 0; kP < nP; kP++)
{
IFeature pF = featureLayer.FeatureClass.GetFeature(kP);
xF[kP] = pF.Extent.XMax; // after help from Martin Maretta
yF[kP] = pF.Extent.YMax;
Debug.Print("Point " + kP.ToString() + " at (" + xF[kP].ToString() + "," + yF[kP].ToString() + ")");
}

// convert coordinates (xF,yF) to raster spatial reference (xR,yR) -----------------------------------
double[] xR = new double[nP];
double[] yR = new double[nP];
IFields pFields = featureLayer.FeatureClass.Fields;
ISpatialReference sr1 = pFields.Field[pFields.FindField(featureLayer.FeatureClass.ShapeFieldName)].GeometryDef.SpatialReference;
ISpatialReference sr2 = rasterProps.SpatialReference;
if (sr1.Name.Equals(sr2.Name,StringComparison.OrdinalIgnoreCase))
{
// spatial references are the same => no need of reprojection
xR = xF;
yR = yF;
}
else // after help from Jamie P. from ESRI Support Team
{
//Point to project
IPoint point = new PointClass() as IPoint;
for (int kP = 0; kP < nP; kP++)
{
point.PutCoords(xF[kP], yF[kP]);
//Geometry Interface to do actual project
IGeometry geometry = point;
geometry.SpatialReference = sr1;
geometry.Project(sr2);
point = geometry as IPoint;
point.QueryCoords(out xR[kP], out yR[kP]);
}
}

// convert from coordinates (xR,yR) to grid cell indexes (kxR,kyR) in raster file --------------------------
IPnt cellSize = rasterProps.MeanCellSize();
kxR = new long[nP];
kyR = new long[nP];
int kI = 0;
for (int kP = 0; kP < nP; kP++)
{
kxR[kI] = (int)Math.Floor((xR[kP] - rasterProps.Extent.XMin) / cellSize.X); // x index starts at XMin
kyR[kI] = (int)Math.Floor((rasterProps.Extent.YMax - yR[kP]) / cellSize.Y); // y index starts at YMax
if (0 <= kxR[kI] && kxR[kI] < rasterProps.Width &&
0 <= kyR[kI] && kyR[kI] < rasterProps.Height)
{
Debug.Print("Point " + kP.ToString() + " indexes= (" + kxR[kI].ToString() + "," + kyR[kI].ToString() + ") => Inside raster area");
kI++; // valid index
}
else
{
Debug.Print("Point " + kP.ToString() + " indexes= (" + kxR[kI].ToString() + "," + kyR[kI].ToString() + ") => Outside raster area");
}
}
Debug.Print("Number of points inside raster area= " + kI.ToString());

if (kI < nP)
{
// some points were discarded => TODO: remove unused array space
nP = kI;
}

return nP;
}
Highlighted
Regular Contributor II
In the VBA examples for 10.0 there was a great example: http://help.arcgis.com/en/sdk/10.0/vba_desktop/conceptualhelp/0001/000100000054000000.htm
You can try to translate it.

Have fun