Select to view content in your preferred language

How to count the adjacent cells that have a different value in a raster dataset?

9505
42
08-08-2018 07:18 PM
MarkLimb
Occasional Contributor

Hi,

I'm using a raster dataset and I'm trying to determine a value for land use adjacency as used by Stathakis & Tsilimigkas (2014) (Measuring the compactness of European medium-sized cities by spatial metrics based on fused data sets. International Journal of Image and Data Fusion, 6(1), 42-64. doi:10.1080/19479832.2014.941018)

This involves counting the number of adjacent cells that have a value different to the centre cell. The attached image illustrates this concept.

It needs to return a new raster dataset where the value of each cell equals this count.

Can anyone suggest a method for doing this in ArcMap 10?

Thank you!

EDIT:

Dan, you're an absolute legend! Amazing dedication and perseverance to assist a total stranger on the other side of the world. Thank you so much for all your effort on this!

Just to update, the main issues experienced had much to do with a combination of differences in arcmap 10.3's and ArcGIS Pro's arcpy, and some major changes in numpy that is used by python 2.7 and python 3.6. Dan worked through them and ended up producing the Toolbox/Script I've attached to this edited post. As a final tip, make sure you run the script in a blank dataframe. You will need to manually set the projection of the final output.

0 Kudos
42 Replies
DanPatterson_Retired
MVP Emeritus

manure... you caught me editing... open in a script editor and delete lines 137 and 138 and save it.

If you get it working.... I am in the process of putting together the whole package (with other stuff... so don't use those yet)

MarkLimb
Occasional Contributor

No worries, mate. I'll wait until you're done.

Many thanks for this... seriously above and beyond!!

0 Kudos
DanPatterson_Retired
MVP Emeritus

Mark... hope I got the right toolbox and script.... 

don't know if the other tools will work, but focus on that focal_stats one

PS... I noticed some funkiness with uploading so that may be part of the problem with things getting stripped out

0 Kudos
MarkLimb
Occasional Contributor

Sorry mate, still an issue:

Executing: FocalStats GSV_Data_Updated_EntClassic_2016 "C:\Research\GIS\Activity Centres.gdb\Adjacency_GSV_Data_Updated_EntClassic_2016" "focal difference"
Start Time: Tue Aug 14 14:29:42 2018
Running script FocalStats...
Failed script FocalStats...

Traceback (most recent call last):
File "C:\Research\GIS\ProximityScript\Attmept4\Scripts\focal_stats.py", line 110, in <module>
desc = arcpy.da.Describe(in_raster)
AttributeError: 'module' object has no attribute 'Describe'

Failed to execute (FocalStats).
Failed at Tue Aug 14 14:29:45 2018 (Elapsed Time: 2.75 seconds)

0 Kudos
DanPatterson_Retired
MVP Emeritus

Mark... due to python 2.6 and old arcpy.  I am using python 3.6.x and the newer arcpy which uses the Data Access module and dictionaries.

Change lines 110-115 to read...

desc = arcpy.Describe(in_raster)
cell_hght = desc.meanCellHeight
cell_wdth = desc.meanCellWidth
cols = desc.width
rows = desc.height
nodata = desc.noDataValue
xtent = desc.extent
pntLL = xtent.lowerLeft
SR = desc.spatialReference
and save the script (you can use Notepad or a pure text editor or if you have a python ide. 
That way, I don't need to roll back the development script
(add familiarity with python programming to your CV )
0 Kudos
MarkLimb
Occasional Contributor

Haha, CV updated!

Unfortunately it returned a number of new errors:

Executing: FocalStats GSV_Data_Updated_EntClassic_2016 "C:\Research\GIS\Activity Centres.gdb\Adjacency_GSV_EntClassic_2016" "focal difference"
Start Time: Tue Aug 14 14:57:31 2018
Running script FocalStats...
Failed script FocalStats...

Traceback (most recent call last):
File "C:\Research\GIS\ProximityScript\Attmept4\Scripts\focal_stats.py", line 125, in <module>
result, nd_out = focal_diff(a, nodata=nodata)
File "C:\Research\GIS\ProximityScript\Attmept4\Scripts\focal_stats.py", line 70, in focal_diff
constant_values=((nodata, nodata), (nodata, nodata)))
File "C:\Python27\ArcGIS10.3\lib\site-packages\numpy\lib\arraypad.py", line 779, in pad
kwargs)
File "C:\Python27\ArcGIS10.3\lib\site-packages\numpy\lib\shape_base.py", line 80, in apply_along_axis
res = func1d(arr[tuple(i.tolist())],*args)
File "C:\Python27\ArcGIS10.3\lib\site-packages\numpy\lib\arraypad.py", line 314, in _constant
return _create_vector(vector, pad_tuple, nconstant[0], nconstant[1])
File "C:\Python27\ArcGIS10.3\lib\site-packages\numpy\lib\arraypad.py", line 37, in _create_vector
vector[:pad_tuple[0]] = before_val
TypeError: long() argument must be a string or a number, not 'NoneType'

Failed to execute (FocalStats).
Failed at Tue Aug 14 14:57:34 2018 (Elapsed Time: 3.20 seconds)

0 Kudos
DanPatterson_Retired
MVP Emeritus

hmmmm never seen a raster without a nodata value....

Try adding....

line 120 and 121...

if nodata is None:
    nodata = -9

just before 

a = arcpy.RasterToNumPyArray….

if that doesn't catch it, I will have to look at the history of np.pad

0 Kudos
MarkLimb
Occasional Contributor

Fewer issues but still a couple...

Executing: FocalStats GSV_Data_Updated_EntClassic_2016 "C:\Research\GIS\Activity Centres.gdb\Adjacency_EntClassic_2016" "focal difference"
Start Time: Tue Aug 14 15:56:41 2018
Running script FocalStats...
Failed script FocalStats...

Traceback (most recent call last):
File "C:\Research\GIS\ProximityScript\Attmept4\Scripts\focal_stats.py", line 127, in <module>
result, nd_out = focal_diff(a, nodata=nodata)
File "C:\Research\GIS\ProximityScript\Attmept4\Scripts\focal_stats.py", line 71, in focal_diff
a_s = stride(ap, win=(3, 3), stepby=(1, 1))
File "C:\Research\GIS\ProximityScript\Attmept4\Scripts\focal_stats.py", line 61, in stride
a_s = as_strided(a, shape=newshape, strides=newstrides, subok=True).squeeze()
TypeError: as_strided() got an unexpected keyword argument 'subok'

Failed to execute (FocalStats).
Failed at Tue Aug 14 15:56:45 2018 (Elapsed Time: 4.23 seconds)

0 Kudos
DanPatterson_Retired
MVP Emeritus

a_s = as_strided(a, shape=newshape, strides=newstride).squeeze()

that one is easy... dump the subok=True... that was a definite change, and a good one... hope it doesn't foul up the next bit

0 Kudos
MarkLimb
Occasional Contributor

Still not there I'm afraid:

Executing: FocalStats GSV_Data_Updated_EntClassic_2016 "C:\Research\GIS\Activity Centres.gdb\Adjacency_EntClassic_2016" "focal difference"
Start Time: Tue Aug 14 16:24:10 2018
Running script FocalStats...
Failed script FocalStats...

Traceback (most recent call last):
File "C:\Research\GIS\ProximityScript\Attmept4\Scripts\focal_stats.py", line 127, in <module>
result, nd_out = focal_diff(a, nodata=nodata)
File "C:\Research\GIS\ProximityScript\Attmept4\Scripts\focal_stats.py", line 71, in focal_diff
a_s = stride(ap, win=(3, 3), stepby=(1, 1))
File "C:\Research\GIS\ProximityScript\Attmept4\Scripts\focal_stats.py", line 61, in stride
a_s = as_strided(a, shape=newshape, strides=newstride).squeeze()
NameError: global name 'newstride' is not defined

Failed to execute (FocalStats).
Failed at Tue Aug 14 16:24:14 2018 (Elapsed Time: 3.34 seconds)

0 Kudos