FRAGSTATS says Reclassify created non-square grids for some of my IMG rasters

410
0
02-21-2014 05:56 PM
MParrish
New Contributor
I'm running ArcGIS 10.2 with an advanced license (student eval).  I have a set of 1188 shapefiles, each containing polygons delineating land cover on a given study site.  I used the PolygonToRaster_conversion tool to convert each shapefile into an IMG raster.  If it matters, each polygon shapefile was split out from a single shapefile that originally contained all of the 1188 sites, and was digitized in ArcMap 9.3.

arcpy.PolygonToRaster_conversion(feature, 'CoverType', 'S:\\perpetual\\fragstats_grids\\ext_bounds_grids\\out_ftr\\' + feature[:-4] + ".img", 'MAXIMUM_AREA', "#", 0.50)


So far so good.  Next, I used Reclassify to change all the various land cover class values to '1', with NODATA for anything outside the external boundary of the study site.  This is because I was interested in doing something with just the shape of the study sites for a question I had. 

outReclass = Reclassify(raster, "Value", RemapRange ([[0,15,1]]), "NODATA")
outReclass.save("S:\\perpetual\\fragstats_grids\\ext_bounds_grids\\out_rascon\\MU_" + raster)


Problem: For 3 of the 1188 rasters, when I import them into FRAGSTATS 4.1, I get an error message saying that FRAGSTATS cannot handle a raster with cells that are non-square.  The other 1185 rasters don't exhibit the same issue.

I tried using Con instead of Reclassify, but got the same error message in FRAGSTATS.  I tried making TIF rasters instead of IMG, but got the same non-square cell error.  I also tried Resample_management to resize the cells, same problem.

outCon = Con(raster, 1, "", "VALUE > 0")
outCon.save("S:\\perpetual\\fragstats_grids\\ext_bounds_grids\\out_rascon\\MU_" + raster)


FRAGSTATS will accept the original (non-reclassified) rasters for those 3 sites.  What I ended up doing was using the PolygonToRaster tool, but using the site ID number as the class, basically skipping the conversion step.  But I'm still bothered by why the non-square cell error appeared.  I assume it's some kind of rounding issue, but if the initial raster was ok, why would Con or Reclassify change the raster cells themselves instead of just reassigning values to each cell?


Here's my entire python code.  I'm new at python, and I may not be doing everything optimally, but I don't think I've made an actual error here.

print ('\n' * 2)
import datetime
import arcpy

"""
ftr()
Reads a directory containing each management unit's land cover shapefile (one shapefile per MU).
Converts all shapefiles into grids (ERDAS IMG format) of a specified cell size.
"""

def ftr():

    # ---------------------------------
    print str(datetime.datetime.now())  + ": Geoprocessing started."
    print ('\n')

    from arcpy import env
    from arcpy.sa import *
    arcpy.CheckOutExtension ('spatial')
    env.workspace = "S:\\perpetual\\fragstats_grids\\ext_bounds_grids\\split_shapefiles_2013_0512"
    env.cellsize = '0.50'

    count = 0
    #featurelist can either look up all features in the workspace, or it can be set to a specific list of features.
    featurelist = arcpy.ListFeatureClasses()
    #featurelist = ['M0072.shp', 'M0098.shp', 'M0422.shp']  
    
    try:
        for feature in featurelist:
            if count < 1189:
                print str(datetime.datetime.now())  + ": " + feature
                arcpy.PolygonToRaster_conversion(feature, 'SiteNum_MP', 'S:\\perpetual\\fragstats_grids\\ext_bounds_grids\\out_ftr\\' + feature[:-4] + ".img", 'MAXIMUM_AREA', "#", 0.50)
                count = count + 1

        print ('\n')
        print str(datetime.datetime.now())  + ": Geoprocessing completed."

    except KeyboardInterrupt:
        print str(datetime.datetime.now())  + ': Run interrupted by user!'
        
    arcpy.CheckInExtension ('spatial')
    


"""
rascon()
Takes rasters containing land cover and converts anything with land cover into a value of 1 and 
everything else as a nodata.  Essentially, makes a solid raster of each MU that can be used in 
fragstats to calculate SHAPE INDEX for each management unit's external bounds.
"""

def rascon():
    print str(datetime.datetime.now())  + ": Geoprocessing started."
    print ('\n')

    from arcpy import env
    from arcpy.sa import *
    arcpy.CheckOutExtension ('spatial')
    env.workspace = "S:\\perpetual\\fragstats_grids\\ext_bounds_grids\\out_ftr"
    env.cellsize = '0.50'
    count = 0
    #raslist can either look up all rasters in the workspace, or it can be set to a specific list of rasters.
    raslist = arcpy.ListRasters()
    #raslist = ['M0072.img', 'M0098.img', 'M0422.img']  

    try:
        for raster in raslist:
            if count < 4:
                print str(datetime.datetime.now())  + ": " + raster
                #outCon = Con(raster, 1.0,)
                #outCon = Con(raster, 1, "", "VALUE > 0")
                #outCon.save("S:\\perpetual\\fragstats_grids\\ext_bounds_grids\\out_rascon\\MU_" + raster)

                outReclass = Reclassify(raster, "Value", RemapRange ([[0,15,1]]), "NODATA")
                outReclass.save("S:\\perpetual\\fragstats_grids\\ext_bounds_grids\\out_rascon\\MU_" + raster)

                count = count + 1

        print ('\n')
        print str(datetime.datetime.now())  + ": Geoprocessing completed."

    except KeyboardInterrupt:
        print str(datetime.datetime.now())  + ': Run interrupted by user!'


    arcpy.CheckInExtension ('spatial')

ftr()
rascon()
    
0 Kudos
0 Replies