topic Re: Zonal statistics - percentiles (reprise) in ArcGIS Spatial Analyst Questions

Zonal statistics - percentiles (reprise)

AndriesPotgieter — Fri, 21 Dec 2018 02:35:11 GMT

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

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

Re: Zonal statistics - percentiles (reprise)

XanderBakker — Fri, 21 Dec 2018 12:04:19 GMT

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-lines-using-arcpy ). Same goes here, without access to the data I will not be able to reproduce the error.

Re: Zonal statistics - percentiles (reprise)

AndriesPotgieter — Fri, 21 Dec 2018 23:32:49 GMT

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

Re: Zonal statistics - percentiles (reprise)

XanderBakker — Sun, 23 Dec 2018 15:52:01 GMT

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.

Re: Zonal statistics - percentiles (reprise)

AndriesPotgieter — Thu, 03 Jan 2019 04:24:56 GMT

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

Re: Zonal statistics - percentiles (reprise)

XanderBakker — Sun, 12 Dec 2021 04:41:59 GMT

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)

Re: Zonal statistics - percentiles (reprise)

AndriesPotgieter — Fri, 04 Jan 2019 05:01:24 GMT

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

Re: Zonal statistics - percentiles (reprise)

XanderBakker — Sat, 05 Jan 2019 21:15:43 GMT

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

Re: Zonal statistics - percentiles (reprise)

AndriesPotgieter — Tue, 08 Jan 2019 00:06:59 GMT

Hi Xander Bakker‌

is it possible to provide the code line for creating 1 single shapefile or table for the resulting zonal statistics here?

thanks

andires

Re: Zonal statistics - percentiles (reprise)

AndriesPotgieter — Tue, 08 Jan 2019 00:09:11 GMT

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()

Re: Zonal statistics - percentiles (reprise)

AndriesPotgieter — Tue, 08 Jan 2019 04:00:26 GMT

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

Re: Zonal statistics - percentiles (reprise)

XanderBakker — Tue, 08 Jan 2019 12:00:23 GMT

Hi Andries Potgieter ,

When you post code in GeoNet, please use these instructions: Posting code with Syntax Highlighting on GeoNet to avoid loosing the indentation which makes the code unreadable.

As far as I understand from the error message there is no Value Attribute Table (VAT) for the raster "c:\PROGRA~2\arcgis\DESKTO~1.5\bin\pol83" and for some reason this file in not readable. Without the data it will be difficult to trace what is going wrong.

Re: Zonal statistics - percentiles (reprise)

XanderBakker — Tue, 08 Jan 2019 12:01:52 GMT

This error seems to be related to the previous one. Please check that the rasters are integer. If they are not, they will not have a VAT which is required in the process. Did you pre-process the rasters the same way as you did with the first one?

Re: Zonal statistics - percentiles (reprise)

XanderBakker — Tue, 08 Jan 2019 12:04:54 GMT

When you want to combine the percentile values for multiple rasters into the same polygon featureclass, the process should be changed drastically, since each raster should yield its own percentile fields each with its own name. At this point it is easier to run the code for each raster, but change the output fields before the next run. This can be done altering the variable fld_prefix for each run.

Re: Zonal statistics - percentiles (reprise)

AndriesPotgieter — Wed, 09 Jan 2019 04:49:20 GMT

Dear Xander Bakker‌

Thanks for the advice. I have changed the code accordingly and it is working now.

Will keep you posted when it has finished.

regards

Re: Zonal statistics - percentiles (reprise)

AndriesPotgieter — Thu, 10 Jan 2019 02:20:05 GMT

Hi Xander Bakker‌,

I have run the Zonal percentile code you wrote on my data for 2 rasters and a snippet of the results is attached as an image here FYI.

thanks again for your help in this.

much appreciated.

Andries

Re: Zonal statistics - percentiles (reprise)

XanderBakker — Thu, 10 Jan 2019 12:48:32 GMT

Glad to hear that it is working! Can you mark the post that answered your question as the correct answer?

Re: Zonal statistics - percentiles (reprise)

AndriesPotgieter — Fri, 01 Mar 2019 03:18:33 GMT

Dear Xander Bakker‌

I am trying to run the python script outside ArcGIS. It takes forever to run in ArcGIS so wanted to speed things up too.

I am getting the following error when executing the first line from command line terminal in windows:

A:\Andries_MainData\data\UAV\boxsite_2016_17\Data>C:\Python27\ArcGIS10.5\python.exe raster_percentile.py
Traceback (most recent call last):
File "raster_percentile.py", line 182, in <module>
main()
File "raster_percentile.py", line 23, in main
arcpy.CreateTable_management(tbl_ws, tbl_name)
File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\management.py", line 17411, in CreateTable
raise e
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Table Location: Dataset A:/data/UAV/boxsite_2016_17/Data/Geodatabase.gdb does not exist or is not supported
Failed to execute (CreateTable).

Can you please advise what I am doing wrong?

I will upload the script also

regards

andries