[GIS] Extract raster by mask output does not completely overlap mask

arcgis-10.0arcpyextract-by-maskrasterspatial-analyst

I have a flow accumulation grid and a shapefile containing catchment areas. Now I want to do some calculations on them so I need to extract data from cells in the raster adjacent to certain catchment areas.

In that process I try to use Extract by Mask on my raster, with one of the catchment polygons as mask (in a loop later on, but now I just need it to work once).

The problem is that the extracted raster is shifted half a cell to the left, compared to both the original raster and the masking polygon. This is messing up the rest of my process, since I try to convert the raster to points and intersect with other polygons later on. And unless the extracted raster is overlapping the original one, I wont get the right polygon.

The original raster and the polygons are aligning perfectly, since the polygons are created from the raster. I have tried to check if there somehow have been some kind of change in the projections, but can't find any. I have also tried in new ArcMap documents, restarting everything several times (this solves the problem surprisingly often) and executing the tools from the python prompt instead of calling my script. I have even tried to execute the Extract by Mask tool from the Toolbox, with same result.

This is the (very simple) code that I use in the python prompt in ArcMap:

inpgs = r'C:\GIS\Catchm_ex.shp'
flowacc = r'C:\GIS\flowacc'
arcpy.MakeFeatureLayer_management(inpgs, 'currentMask', '"HydroID" = 269634')
outRaster = arcpy.sa.ExtractByMask(flowacc, 'currentMask')

I have some sample files here

I use Windows 7, ArcInfo 10.0 and Python 2.6.

I have tried Extract by Polygon now, and that seems to work fine. Is it just my machine malfunctioning on ExtractByMask?
Only drawback is that I'm having problems using that function in my script, since I can't work out how to get the coordinates of my masking polygon (only covering 9 cells, with a point in center if that helps).

I have now got it working with Extract by polygon by using the code below (modified from Esri's support pages). Unfortunately, this does not solve my problem. It appears as if the output (from this tool as well as Extract by mask) is sometimes shifted half a cell to the left, and sometimes not. Mostly it is, but occasionally it actually is located where it should be.

So now I'm really confused. Is this just happening to me, on my computer? If anyone else have experienced the same, please write a comment so I know that I'm not all alone 🙂 To me this problem seems to be a too big of an issue to be a bug, which leaves me with this being a problem that I have created myself along the way.

Any ideas on how to avoid it?

def ReadPoints(infc):
    try:
        import arcpy, numpy
    except ImportError:
        sys.exit('Unable to import arcpy')

    # Identify the geometry field
    desc = arcpy.Describe(infc)
    shapefieldname = desc.ShapeFieldName
    points = []

    # Create search cursor
    rows = arcpy.SearchCursor(infc)

    # Enter for loop for each feature/row
    for row in rows:
        # Create the geometry object
        feat = row.getValue(shapefieldname)
        partnum = 0

        # Step through each part of the feature
        for part in feat:
            # Step through each vertex in the feature
            for pnt in feat.getPart(partnum):
                if pnt:
                    points.append([int(pnt.X), int(pnt.Y)])

    pt = numpy.array(points)
    pts = []
    for i in range(len(pt)):
        t = pt[i]
        pts.append(arcpy.Point(t[0],t[1]))
    return pts

Best Answer

This may not be helpful now, but this is a known issue with the Extract by Mask tool that was obviously not fixed for ArcGIS 10: http://support.esri.com/en/knowledgebase/techarticles/detail/34448

I don't know the scripting solution . . .

Related Question