0

I need to reclassify two raster sets and turn them into polygons for further processing. I have to conduct multiple classifications so doing it manually is not a real option. Also I am bound to doing it in ArcGIS/ArcPy.

However, if I reclass raster1 (Paris) with ArcPy it returns an empty raster (no Error Codes or anything) and of course the conversion to polygons fails right after. When doing it in ArcGIS with some random classes it works perfectly fine. Also raster 2 (Graz) works fine with my code. I assume it's a problem with the data, since raster 1 contains some NAs. But the same problem occurs, if I replace NAs with 0.

Did any of you ever run in similar problems?

If you want to try it on your own, I uploaded a snippet of the data here and attached a code snippet below.

The remap list currently looks like this, but will vary with the number of classes. Since this format works fine for Raster 2, I don't think that this is the problem, but feel free to prove me wrong

For Raster 1:

[[0, 1, 'NODATA'], [1, 116, 1], [232, 348, 2], [348, 464, 3], [464, 580, 4], [580, 696, 5], [696, 812, 6], [812, 928, 7], [928, 1044, 8], [1044, 1160, 9]]

For Raster 2:

[[329, 359, 0], [359, 389, 1], [389, 419, 2], [419, 449, 3], [449, 479, 4], [479, 509, 5], [509, 539, 6], [539, 569, 7], [569, 599, 8], [599, 629, 9], [629, 659, 10]]

Even though I have 5+ years of experience, I can't determine the problem here.

import arcpy as arc
import math as ma


# set environment settings
arc.CheckOutExtension('Spatial')


def getMinMax(inRaster):

    mini = arc.GetRasterProperties_management(inRaster, "MINIMUM")
    maxi = arc.GetRasterProperties_management(inRaster, "MAXIMUM")
    mini = ma.floor(float(mini.getOutput(0).replace(',', '.')))
    maxi = ma.ceil(float(maxi.getOutput(0).replace(',', '.')))

    return mini, maxi

def createRemapList(minVal, maxVal, stepSize):
    valueList = range(int(minVal), int(maxVal+stepSize), int(stepSize))
    mapList = []
    for count, val in enumerate(valueList):
        if count == len(valueList) - 1:
            break
        if val == 0:
            listAttach = [0, 1, 'NODATA']
        elif valueList[count-1] == 0:
            listAttach = [1, valueList[count], count]
        else:
            listAttach = [valueList[count], valueList[count+1], count]

        mapList.append(listAttach)
    print(mapList)

    return mapList

pop_path = './Data/Pop/Paris_con_reclass.tif'
dgm_path = './Data/DGM/DGM_Graz_sub_sub.tif'


# reclassifying rasters to stratify data into groups
# get min and max value of data
minPop, maxPop = getMinMax(pop_path)
minDGM, maxDGM = getMinMax(dgm_path)
# define number of classes
numClasses = 10
# define equal steps between values
StepsPop = round((maxPop - minPop)/numClasses)
StepsDGM = round((maxDGM - minDGM)/numClasses)

reclassPathPop = "./Data/Pop/Reclass_Paris.tif"
reclassPathDGM = "./Data/DGM/Reclass_Graz_Sub.tif"

mapListPop = createRemapList(minVal=minPop, maxVal=maxPop, stepSize=StepsPop)
mapListDGM = createRemapList(minVal=minDGM, maxVal=maxDGM, stepSize=StepsDGM)

rasterPop = arc.Raster(pop_path)
rasterDGM = arc.Raster(dgm_path)

pop_reclass = arc.sa.Reclassify(pop_path, "Value",
                                remap=arc.sa.RemapRange(mapListPop))
pop_reclass.save(reclassPathPop)

dgm_reclass = arc.sa.Reclassify(rasterDGM, "Value",
                                remap=arc.sa.RemapRange(mapListDGM))
dgm_reclass.save(reclassPathDGM)

# Transform reclass to polygon
outPolygonReclassDGM = './Data/DGM/Reclass_Graz_Sub_Poly.shp'
arc.RasterToPolygon_conversion(in_raster=reclassPathDGM,
                               out_polygon_features=outPolygonReclassDGM,
                               simplify='NO_SIMPLIFY',
                               create_multipart_features='MULTIPLE_OUTER_PART')

outPolygonReclassPop = './Data/DGM/Reclass_Paris_Poly.shp'
arc.RasterToPolygon_conversion(in_raster=reclassPathPop,
                               out_polygon_features=outPolygonReclassPop,
                               simplify='NO_SIMPLIFY',
                               create_multipart_features='MULTIPLE_OUTER_PART')
2
  • It may or may not matter but are you using ArcPy and Spatial Analyst with ArcGIS Pro or ArcMap?
    – PolyGeo
    Commented May 28, 2021 at 2:37
  • I'm using it with ArcMap. But since this is a very basic tool, I don't think it matters a lot...
    – Thom Spi
    Commented May 28, 2021 at 15:42

1 Answer 1

0

Apparently ArcPy doesn't like to process two different rasters with different extent etc one after the other in the same script... There was something blocked in the environment, when the first raster was read, what made the result of raster 2's reclass empty.

However, splitting processing of both raster into two scripts made it work.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.