Zonal statistics - percentiles (reprise)

2347
17
12-20-2018 06:35 PM
AndriesPotgieter
New Contributor II

Brached from https://community.esri.com/thread/95403 

Hi Xander Bakker 

I have also run your multiple raster percentile code. It runs upto loading the created polygon file:

execfile(r'F:\\data\\UAV\\Lodging_model\\GeoDatabase\\raster_percentile.py')
Processing Raster 'after_raster'
- Processing polygon: 0
- ExtractByPolygon...
- saved raster as ras0pol0
- fill dict with Value x Count
Runtime error
Traceback (most recent call last):
File "<string>", line 1, in <module>
File "F:\\data\\UAV\\Lodging_model\\GeoDatabase\\raster_percentile.py", line 171, in <module>
main()
File "F:\\data\\UAV\\Lodging_model\\GeoDatabase\\raster_percentile.py", line 67, in main
dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(outname, flds_ras)}
RuntimeError: cannot open 'ras0pol0'

here is the first part of the code:

import arcpy
import os

def main():
# settings
ras_ws = r"F:/data/UAV/Lodging_model/GeoDatabase/raster/" # input raster workspace
fc = r"F:/data/UAV/Lodging_model/GeoDatabase/Hermitage_Lodging_Grid.shp" # input polygons
tbl = r"F:/data/UAV/Lodging_model/GeoDatabase/raster/Zonal_Perc_A" # output related table

lst_perc = [2, 5, 10, 90, 95, 98] # list of percentiles to be calculated
fld_prefix = "Perc_"
fld_poloid = "POL_OID"
fld_rasname = "RasterName"

---------------------

not sure why it cant load the ras0pol0 ?  It opens in acrgis no p[roblem from the geodatabase...

cheers

adnries

0 Kudos
17 Replies
XanderBakker
Esri Esteemed Contributor

Hi Andries Potgieter ,

It there a way you can share the data with me? It will be very difficult to discover why there is an error without access to the data. You could also try and use a File Geodatabase (that would be the first thing I would do). 

I also noted that you commented that all percentiles are the same in an unrelated blog here (https://community.esri.com/people/xander_bakker/blog/2014/07/04/japanese-visual-multiplication-with-... ). Same goes here, without access to the data I will not be able to reproduce the error.

0 Kudos
AndriesPotgieter
New Contributor II

Hi Zander

Thanks for your reply. I will upload the two files during the next two days or can i send you a dropbox link?

Andries

0 Kudos
XanderBakker
Esri Esteemed Contributor

Hi a.potgieter@uq.edu.au , 

If you upload the files (<50MB) at this thread other users can chime in and help you too. If the link to dropbox is available for everyone and the files are bigger than 50MB Dropbox is good alternative. 

I must say that I'm a bit busy these days with delivering some products before the end of this year, but if I find a moment I will certainly give it a shot. 

0 Kudos
AndriesPotgieter
New Contributor II

Hi Xander Bakker‌,

Apologies for the late reply.

The files are too large to upload here, unfortunately.

here is the link that has the plot grid data shapefile and the raster. It would be good if you can make it run on more than 1 rater by maybe duplicating the raster.

the files are copied to this dropbox link:

Dropbox - data_Xander - Simplify your life 

I hope you can find the problem.

cheers

andries

0 Kudos
XanderBakker
Esri Esteemed Contributor

Hi Andries Potgieter ,

I only have access to ArcGIS Pro here, so this is what I did:

  • I converted the shapefile to a featureclass in a File Geodatabase and assigned the same projection as the raster (since they did not match)
  • The raster is a floating raster (this will only work with integer raster), so I use the math algebra to create and integer raster multiplying the input with 100.
  • I changed the script a bit and included an additional Int statement to be sure that every extracted raster was going to be integer as well (see script below)
  • I loaded the integer raster and the featureclass into ArcGIS Pro and ran this script:

import arcpy  
  
def main():  
    # settings  
    ras = "raster_int2.tif" # r"C:\Users\xande\OneDrive\Documenten\ArcGIS\Projects\PercentileRaster\PercentileRaster.gdb\raster_int2"  
    fc = "lodging_grid" # r"C:\Users\xande\OneDrive\Documenten\ArcGIS\Projects\PercentileRaster\PercentileRaster.gdb\lodging_grid"  
    lst_perc = [2, 5, 10, 25, 50, 75, 90, 95, 98] # list of percentiles to be calculated  
    fld_prefix = "Perc_"  
    start_at = 0  
  
    arcpy.env.workspace = r"C:\Users\xande\OneDrive\Documenten\ArcGIS\Projects\PercentileRaster\PercentileRaster.gdb" # "IN_MEMORY"  
    arcpy.env.overwriteOutput = True  
  
    # add fields if they do not exist in fc  
    # be aware, fields that exist will be overwritten!  
    print("filling field list and adding fields")
    flds = []  
    for perc in lst_perc:  
        fld_perc = "{0}{1}".format(fld_prefix, perc)  
        flds.append(fld_perc)  
        if not FieldExist(fc, fld_perc):  
            arcpy.AddField_management(fc, fld_perc, "LONG")  
    print("flds={0}".format(flds))
  
    # Enable Spatial analyst  
    arcpy.CheckOutExtension("Spatial")  
  
    # get cell area  
    desc = arcpy.Describe(ras)  
    cell_area = desc.meanCellHeight * desc.meanCellWidth  
  
    # loop through polygons  
    print("loop through polygons")
    flds.append("SHAPE@")  
    i = 0  
    with arcpy.da.UpdateCursor(fc, flds) as curs:  
        for row in curs:  
            i += 1  
            print("Processing polygon: {0}".format(i))
            if i >= start_at:  
                polygon = row[flds.index("SHAPE@")]  
  
                if polygon.area < cell_area:  
                    polygon = polygon.buffer((desc.meanCellHeight + desc.meanCellWidth) / 4)  
  
                # Execute ExtractByPolygon (you can't send the polygon object)  
                print(" - ExtractByMask...")
  
                # arcpy.sa.ExtractByPolygon(ras, lst_parts, "INSIDE")  
                ras_pol1 = arcpy.sa.ExtractByMask(ras, polygon)  
                del polygon  
                outname = "rass{0}".format(i)  
                ras_pol1.save(outname)  
                print(" - saved raster1 as {0}".format(outname))
                    
                print("make integer...")
                ras_pol = arcpy.sa.Int(ras_pol1)
                outname = "ras{0}".format(i)  
                ras_pol.save(outname)  
                print(" - saved raster as {0}".format(outname))
  
                # create dictionary with value vs count  
                print(" - fill dict with Value x Count")  
                flds_ras = ("VALUE", "COUNT")  
                dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(outname, flds_ras)}  
                # del ras_pol  
  
                # calculate number of pixels in raster  
                print(" - determine sum")
                cnt_sum = sum(dct.values())  
                print(" - sum={0}".format(cnt_sum))
  
                # loop through dictionary and create new dictionary with val vs percentile  
                print(" - create percentile dict")
                dct_per = {}  
                cnt_i = 0  
                for val in sorted(dct.keys()):  
                    cnt_i += dct[val]  
                    dct_per[val] = cnt_i / cnt_sum  
                del dct  
  
                # loop through list of percentiles  
                print(" - iterate perceniles")
                for perc in lst_perc:  
                    # use dct_per to determine percentiles  
                    perc_dec = float(perc) / 100  
                    print("  - Perc_dec for is {0}".format(perc_dec))
                    pixval =  GetPixelValueForPercentile(dct_per, perc_dec)  
                    print("  - Perc for {0}% is {1}".format(perc, pixval))
  
                    # write pixel value to percentile field  
                    print("  - Store value")
                    fld_perc = "{0}{1}".format(fld_prefix, perc)  
                    row[flds.index(fld_perc)] = pixval  
                del dct_per  
  
                # update row  
                print(" - update row")
                curs.updateRow(row)  
  
    # return SA license  
    arcpy.CheckInExtension("Spatial")  
  
def GetPixelValueForPercentile(dctper, percentile):  
    """will return last pixel value 
      where percentile LE searched percentile."""  
    pix_val = -1  
    try:  
        pix_val = sorted(dctper.keys())[0] # initially assign lowest pixel value  
        for k in sorted(dctper.keys()):  
            perc = dctper[k]  
            if perc <= percentile:  
                pix_val = k  
            else:  
                break  
    except Exception as e:  
        pass  
    return pix_val  
  
def FieldExist(featureclass, fieldname):  
    """Check if field exists"""  
    fieldList = arcpy.ListFields(featureclass, fieldname)  
    return len(fieldList) == 1  
  
if __name__ == '__main__':  
    main()  ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

This seemed to work. I will attach the resulting featureclass (as shapefile)

0 Kudos
AndriesPotgieter
New Contributor II

 Magic!! Xander Bakker

I just multiplied the raster height by 1000 at your INt convert line

ras_pol = arcpy.sa.Int(ras_pol1 * 1000)

I would like to run it on 2 rasters and have all data in a table format for each plot. I assume I can copy the second raster to geodatabase and include the INT code transform line in it?

much appreciated

Andries

0 Kudos
XanderBakker
Esri Esteemed Contributor

Yes, that is certainly possible. Great to hear that it is working!

0 Kudos
AndriesPotgieter
New Contributor II

Hi Xander Bakker

I have adjusted the paths fro using your code of creating zonal percentiles on 2 rasters or more:

I am getting the following error I suspect it is a matter of wrong path somewhere/

can you please advise?

regards

andries

error:

Runtime error Traceback (most recent call last): File "<string>", line 1, in <module> File "F:/data/UAV/Lodging_model/Geodatabase/Raster_percentile.py", line 171, in <module> main() File "F:/data/UAV/Lodging_model/Geodatabase/Raster_percentile.py", line 59, in main ras_pol = arcpy.sa.ExtractByPolygon(ras, lst_parts, "INSIDE") File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Functions.py", line 1482, in ExtractByPolygon extraction_area) File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Utils.py", line 53, in swapper result = wrapper(*args, **kwargs) File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Functions.py", line 1477, in Wrapper extraction_area) File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\geoprocessing\_base.py", line 510, in <lambda> return lambda *args: val(*gp_fixargs(args, True)) ExecuteError: ERROR 999999: Error executing function. The table was not found. [VAT_after_raster] ERROR 010060: Could not open and write to c:\PROGRA~2\arcgis\DESKTO~1.5\bin\pol83. Failed to execute (ExtractByPolygon).

code:

import arcpy
import os

def main():
# settings
ras_ws = r"F:/data/UAV/Lodging_model/Percentile.gdb/" # input raster workspace
fc = r"F:/data/UAV/Lodging_model/Percentile.gdb/Hermitage_Lodging_Grid" # input polygons
tbl = r"F:/data/UAV/Lodging_model/Percentile.gdb/Zonal_Perc" # output related table

lst_perc = [2, 98] # list of percentiles to be calculated
fld_prefix = "Perc_"
fld_poloid = "POL_OID"
fld_rasname = "RasterName"

# create table
tbl_ws, tbl_name = os.path.split(tbl)
arcpy.CreateTable_management(tbl_ws, tbl_name)

# add fields and fill fields list
flds_tbl = [fld_poloid, fld_rasname]
arcpy.AddField_management(tbl, fld_poloid, "LONG")
arcpy.AddField_management(tbl, fld_rasname, "TEXT", 50)
for perc in lst_perc:
fld_perc = "{0}{1}".format(fld_prefix, perc)
arcpy.AddField_management(tbl, fld_perc, "LONG")
flds_tbl.append(fld_perc)

# Enable Spatial analyst
arcpy.CheckOutExtension("Spatial")

# get list of rasters
lst_ras = getListOfRasterFromWS(ras_ws)

# environments
arcpy.env.workspace = "F:/data/UAV/Lodging_model/Percentile.gdb/"
arcpy.env.overwriteOutput = True

# create dictionary with polygons OID vs lst_parts
dct_pols = createDictPolygons(fc, fld_poloid)

# start insert cursor for table
with arcpy.da.InsertCursor(tbl, flds_tbl) as curs_tbl:

# start loop through rasters
i = 0
for ras_name in lst_ras:
ras = os.path.join(ras_ws, ras_name)
i += 0
print "Processing Raster '{0}'".format(ras_name)

# loop through polygons
for oid, lst_parts in dct_pols.items():
lst_row = [oid, ras_name]
print " - Processing polygon: {0}".format(oid)

# Execute ExtractByPolygon (you can't send the polygon object)
print " - ExtractByPolygon..."
ras_pol = arcpy.sa.ExtractByPolygon(ras, lst_parts, "INSIDE")
outname = "ras{0}pol{1}".format(i, oid)
ras_pol.save(outname)
print " - saved raster as {0}".format(outname)

# create dictionary with value vs count
print " - fill dict with Value x Count"
flds_ras = ("VALUE", "COUNT")
dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(outname, flds_ras)}

# calculate number of pixels in raster
print " - determine sum"
cnt_sum = sum(dct.values())
print " - sum={0}".format(cnt_sum)

# loop through dictionary and create new dictionary with val vs percentile
print " - create percentile dict"
dct_per = {}
cnt_i = 0
for val in sorted(dct.keys()):
cnt_i += dct[val]
dct_per[val] = cnt_i / cnt_sum

# loop through list of percentiles
print " - iterate percentiles"
for perc in lst_perc:
# use dct_per to determine percentiles
perc_dec = float(perc) / 100
print " - Perc_dec for is {0}".format(perc_dec)
pixval = GetPixelValueForPercentile(dct_per, perc_dec)
print " - Perc for {0}% is {1}".format(perc, pixval)

# write pixel value to percentile field
print " - Store value in list"
# fld_perc = "{0}{1}".format(fld_prefix, perc)
# row[flds.index(fld_perc)] = pixval
lst_row.append(pixval)

# update row
print " - insert row"
row = tuple(lst_row)
curs_tbl.insertRow(row)

# return SA license
arcpy.CheckInExtension("Spatial")
print "Ready..."


def createDictPolygons(fc, fld_poloid):
flds = ("OID@", "SHAPE@", fld_poloid)

# add pol oid field for link
if not FieldExist(fc, fld_poloid):
arcpy.AddField_management(fc, fld_poloid, "LONG")

dct_pols = {}
with arcpy.da.UpdateCursor(fc, flds) as curs:
for row in curs:
oid = row[0]
polygon = row[1]
row[2] = oid

lst_parts = []
if polygon.partCount == 1:
for part in polygon:
for pnt in part:
x, y = pnt.X, pnt.Y
lst_parts.append(arcpy.Point(x, y))
else:
for part in polygon:
lst_crds = []
for pnt in part:
x, y = pnt.X, pnt.Y
lst_crds.append(arcpy.Point(x, y))
lst_parts.append(lst_crds)

# add to dict and store pol oid
dct_pols[oid] = lst_parts
curs.updateRow(row)
return dct_pols


def getListOfRasterFromWS(ras_ws):
arcpy.env.workspace = ras_ws
return arcpy.ListRasters()


def GetPixelValueForPercentile(dctper, percentile):
"""will return last pixel value
where percentile LE searched percentile."""
try:
srt_keys = sorted(dctper.keys())
pix_val = srt_keys[0]
for k in srt_keys:
perc = dctper
if perc <= percentile:
pix_val = k
else:
break
return pix_val
except Exception as e:
print " - GetPixelValueForPercentile error: {0}".format(e)
print "dctper:\n{0}".format(dctper)
return -9999

def FieldExist(featureclass, fieldname):
"""Check if field exists"""
import arcpy
fieldList = arcpy.ListFields(featureclass, fieldname)
return len(fieldList) == 1

if __name__ == '__main__':
main()

0 Kudos
AndriesPotgieter
New Contributor II

Hi Xander Bakker

after creating a new geodatabase and copying all relevant rasters and shapefiels to it I now get this error:

execfile(r'F:/data/UAV/Lodging_model/Geodatabase/Raster_percentile.py')
Processing Raster 'after_raster'
- Processing polygon: 1
- ExtractByPolygon...
- saved raster as ras0pol1
- fill dict with Value x Count
Runtime error
Traceback (most recent call last):
File "<string>", line 1, in <module>
File "F:/data/UAV/Lodging_model/Geodatabase/Raster_percentile.py", line 171, in <module>
main()
File "F:/data/UAV/Lodging_model/Geodatabase/Raster_percentile.py", line 67, in main
dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(outname, flds_ras)}
RuntimeError: cannot open 'ras0pol1'
>>>

Would be great if you can perhaps advise me what the likely problem is?

cheers

andries

0 Kudos