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