Count the percentage of positive and Negative pixel within the zone.
I noticed something strange in the data that you attached:
The folder Count_Pixel.gdb is a folder and not a File Geodatabase as I suspected. The folder contains two subfolders, one with a TIF for the NPP slope values and another with an Esri Grid file with the rasterized zones.
To avoid confusion I placed them in a separate folder and adjusted the code to run with the data. See code below:
#-------------------------------------------------------------------------------
# Name: count_pos_neg_pixels.py
# Purpose:
#
# Author: Xander
#
# Created: 06-08-2017
#-------------------------------------------------------------------------------
def main():
print "load modules..."
import arcpy
import numpy as np
arcpy.env.overwriteOutput = True
values_ras = r'C:\GeoNet\CountPixels\NPP_slope.TIF'
zones_ras = r'C:\GeoNet\CountPixels\acz_raster'
nodata = -3.4028235e+38
print "create numpy rasters..."
np_zones = arcpy.RasterToNumPyArray(zones_ras)
np_values = arcpy.RasterToNumPyArray(values_ras)
print "read props numpy raster..."
rows = np_zones.shape[0]
cols = np_zones.shape[1]
print " - rows:", rows
print " - cols:", cols
print "loop through pixels..."
pix_cnt = 0
dct = {}
for row in range(rows):
for col in range(cols):
pix_cnt += 1
if pix_cnt % 5000 == 0:
print " - row:", row, " col:", col, " pixel:", pix_cnt
zone = np_zones[row, col]
value = np_values[row, col]
if zone > 1:
if zone in dct:
stats = dct[zone]
stats = UpdateStats(stats, value, nodata)
else:
stats = CreateStats(value, nodata)
dct[zone] = stats
print "\nZone\tsum\tcnt\tmean\tmin\tmax\tcnt_neg\tcnt_pos\tcnt_nodata"
for zone, stats in sorted(dct.items()):
print zone, '\t', '\t'.join(str(a).replace('.', ',') for a in stats)
def UpdateStats(stats, value, nodata):
# sum, cnt, mean, min, max, cnt_neg, cnt_pos, cnt_nodata
if value > nodata / 2.0:
stats[0] += value
stats[1] += 1
stats[2] = float(stats[0]) / float(stats[1])
if value < stats[3]:
stats[3] = value
if value > stats[4]:
stats[4] = value
if value < 0:
stats[5] += 1
if value > 0:
stats[6] += 1
else:
stats[7] += 1
return stats
def CreateStats(value, nodata):
# sum, cnt, mean, min, max, cnt_neg, cnt_pos, cnt_nodata
if value > nodata / 2.0:
return [value, 1, value, value, value, 1 if value < 0 else 0, 1 if value > 0 else 0, 0]
else:
return [0, 0, None, 99999999, -9999999, 0, 0, 1]
if __name__ == '__main__':
main()The result was loaded into Excel:
When you do it the way Dan and Jayanta explained, you will create a raster with only the positive values:
Create a raster with only the negative values:
Do a Zonal Statistics as Table on the raster with the positive values:
.... and another one on the negative values:
Which will produce this result for the positive values:
and this one for the negative values:
When I compare the count for both methods, they correspond.
Hi Xander Bakker Yes, now it's appearing perfectly and got exact correct pixel calculation. thanks for such wonderful script and it will save the time and also reduce the chances of human introduce error.
And the last column (acz_raster) did you type it manually into excel, based on the zone sequence in the zone raster or it's automatically come through script?
Because for me it's not appearing. it looks like (below attached) after paste into excel, after "cnt_nodata", there is no other column.
The table which you attached there is complete information. but when i am running not getting the same, is there anything need to change?
Thank you
regards
Shouvik Jha
Guilty as charged... I did copied and pasted these values in the last field, from the attribute table of the raster. However, see below a small snippet of code to extract the value into a dictionary. You will need the arcpy.SearchCursor since the arcpy.da.SearchCursor does not support attribute tables of rasters datasets.
For instance my example raster has the following fields:
You can use the following code to create a dictionary:
def main():
import arcpy
zones_ras = r'C:\GeoNet\CountPosNegPixels\data2.gdb\test'
fld_value = "Value"
fld_zone = "Zone"
dct = {}
curs = arcpy.SearchCursor(zones_ras)
for row in curs:
value = row.getValue(fld_value)
zone = row.getValue(fld_zone)
dct[value] = zone
for val, zone in dct.items():
print "value {0} corresponds to {1}".format(val, zone)
if __name__ == '__main__':
main()This will print:
value 1 corresponds to Zone 8
value 2 corresponds to Zone 7
value 3 corresponds to Zone 6
value 4 corresponds to Zone 5
value 5 corresponds to Zone 4
value 6 corresponds to Zone 3
value 7 corresponds to Zone 2
value 8 corresponds to Zone 1In the code after creating the dictionary, when you write the result you can convert the raster value back to zone by using dct[value] to get the zone description.
