[GIS] Calculating first quartile of a raster by zonal statistics

arcpyrasterraster-calculatorspatial-analystzonal statistics

I want to calculate the 1st quartile (25th percentile) of a raster (beneath polygons) using zonal statistics but it does not this option.

Zonal statistics has: MIN, MAX, MEAN, MAJORITY, RANGE, and does not have 1st quartile. How can I do that?

Best Answer

Dear dof1985 thanks for your comment I did it by the following code : 

import arcpy

from arcpy import env

from arcpy.sa import *

import arcpy, os

import arcpy.da

env.workspace = "E:/spring/cases/main/9_jan_2016/New Folder/output"
output1=r"E:\main\9_jan_2016\New Folder"
output=r"E:\main\9_jan_2016\New Folder\output\test"
output2=r"E:\main\9_jan_2016\New Folder\output\reza1.gdb"
output3=r"E:\main\9_jan_2016\New Folder\output"

fc=output3+os.sep+"1800.shp"
inraster=output1+os.sep+"1800TIR1BTC.tif"
fc2=output2+os.sep+"percent07"
g=0
arcpy.CreateTable_management (output2,"percent07",       output3+os.sep+str(1800)+".shp")
arcpy.CheckOutExtension("Spatial")
inraster=Int(inraster)

for row in arcpy.da.SearchCursor(fc, ["FID","SHAPE@"]):
    print(row[0])
    average=[]
    outfc = os.path.join(output, "fc" + str(g))
    arcpy.CopyFeatures_management(row[1], outfc)                    
    # Check out the ArcGIS Spatial Analyst extension license
    arcpy.CheckOutExtension("Spatial")


   # Execute ExtractByPolygon

    TabulateArea(inraster,"Value",output+os.sep+ "fc" + str(g)+".shp","Id",output+os.sep+"reza"+str(g)+".dbf")        
    A=int(arcpy.GetCount_management(output+os.sep+"reza"+str(g)+".dbf").getOutput(0))
    print("A"+str(A))
    counti=A*0.25
    print("counti"+str(counti))
    b=0
    i=0
    if (A!=0):

       arcpy.Sort_management(output+os.sep+"reza"+str(g)+".dbf", output+os.sep+"rezasorted"+str(g)+".dbf",[["Value","ASCENDING"]])
       for row in arcpy.da.SearchCursor(output+os.sep+"rezasorted"+str(g)+".dbf", ["Value"]):
           if  i<=counti:
           #cnt_sum = sum(extPolygonOut.values())
            print("row"+str(row[0]))
            b=b+row[0]
            i=i+1            
           else :
             print("loop  finished")                                           
             break

       a=b/i    


       average.append(a)
       with arcpy.da.InsertCursor(fc2,["TIR1"])as cursor :
                for row in average:

                    cursor.insertRow(average)

       arcpy.Delete_management(output+os.sep+ "fc" + str(g)+".shp")
       arcpy.Delete_management(output+os.sep+"reza"+str(g)+".dbf")
       arcpy.Delete_management(output+os.sep+"rezasorted"+str(g)+".dbf")
       g=g+1
       print(average) 
       #del(average)
    else :
          arcpy.Delete_management(output+os.sep+ "fc" + str(g)+".shp")
          arcpy.Delete_management(output+os.sep+"reza"+str(g)+".dbf")
          arcpy.Delete_management(output+os.sep+"rezasorted"+str(g)+".dbf")
          g=g+1
Related Question