0

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?

  • I guess you are using Arcgis? do you have access to other software? – dof1985 Jul 17 '16 at 15:15
  • Also, see this thread from Geonet; if it works you can post the code as an answer. – dof1985 Jul 17 '16 at 15:19
  • thanks for your consideration. I am using ArcGIS. I looked at that script it is used for percentile not calculation the 1th quarter. I have written some codes using arcpy and it is preferable to use ARCGIS. – Khandan Jul 17 '16 at 18:44
  • 3
    1st quarter = 25th percentile. Change line 25 from lst_perc = [2, 5, 10, 90, 95, 98] to lst_perc = [25] – user2856 Jul 17 '16 at 19:30
  • thanks for your precious comment. Is it possible to calculate the average of 1th quarter of values? – Khandan Jul 17 '16 at 20:51
1

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   
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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