How to find: 1) Pixels in a raster which distance between them is maximized; and 2) have also the maximum value from all the pixels

Question

I am not sure if this is possible, at least I could not find the solution. I have a raster and each pixel has a value. What I try to do is to write a function that finds:

1. The 5 pixels the distance between those is maximized. Say the 5 pixels most far away from each other.
2. And, at the same time, the value of those pixels is the maximum possible while satisfying the first condition.

The raster looks like below, yet this question is general and may be applied to any sort of rasters. Particularly this raster has values on solar radiation (W/m2). http://www.mediafire.com/download/t7xl9adxz3764f3/solar_dni_20130102_00UT.txt

setwd("/Users/mercel/Documents/Australia")
library(SDMTools)
library(raster)
file.rename("solar_dni_20130102_11UT.txt",
        "solar_dni_20130102_11UT.asc")
dat <- read.asc("solar_dni_20130102_11UT.asc")
str(dat)
r <- raster.from.asc(dat)
...

The version I am using is ArcMap 10.2. And I have Python 27 and Python 35 installed. I have tried to do it applying 'Euclidean Distance' and other functions per default in Desktop, yet I believe this may required some sort of optimization coding, which I do not have previous experience with. Hence I cannot provide an example of the code developed so far, because really I don't know even where to start from. Therefore someone with more experience than I, may know how to approach this problem.

Error found by executing the code of FelixIP below, in the Python window of ArcMap 10.2:

Runtime error 
Traceback (most recent call last):
  File "<string>", line 24, in <module>
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\__init__.py", line 202, in SetProgressor
    return gp.setProgressor(type, message, min_range, max_range, step_value)
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\geoprocessing\_base.py", line 83, in setProgressor
    self._gp.SetProgressor(*gp_fixargs(args, True)))
AttributeError: Object: Error in parsing arguments for SetProgressor

Error when copying the .py into C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin, and try to run it from PythonWin 2.7.5 (32):

File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 323, in RunScript
debugger.run(codeObject, __main__.__dict__, start_stepping=0)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\__init__.py", line 60, in run
_GetCurrentDebugger().run(cmd, globals,locals, start_stepping)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\debugger.py", line 654, in run
exec cmd in globals, locals
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\Script.py", line 1, in <module>
try: arcpy
NameError: name 'time' is not defined`

Answer 1

I took this raster, defined it’s projection as GCS_WGS_1984 and projected it to GDA_1994_Geoscience_Australia_Lambert. It might not be the best projection to use over such vast extent. I think you should post separate question re best projection to compare distances all over the continent.

I re-sampled it to 125 km*125 km cell size and converted results to point shapefile containing 571 points total. I placed this shapefile at the top of my table of content in mxd.

I defined total threshold for total as (mean + 0.5*standard deviations)*n, where mean and standard deviation are statistics of field “GRID_CODE” in above points and n is number of points to be found. See relevant lines in script below.

I ran script from active mxd as shown, i.e. n=3.

RESULT:

Warning! It took close to 40(!) minutes to complete task for 3(!) points. So use at your own risk if you’d like to pick 5 points, but first have a look at http://stattrek.com/online-calculator/combinations-permutations.aspx

It will give you an idea regarding count of combinations to go through to do the job, because this is exactly what script does. It calculates totals for n points and rejects them, if it is below threshold. Otherwise it calculates the distance between pairs of points and proceeds with next combination of points. Finally, it selects n points with largest distance between them.

Note that very low threshold will give you points on the perimeter (similar to result shown), while very high one will result in points crowded around ‘high point’ or script failure.

On the bright side of things:

• I used my very old PC at home to test it
• Script has progress bar, so click Cancel at any time if you see that it’s going to take few human lives to complete.

To locate points more accurately, you can always find highest points in original raster within 125 km of points selected.

SCRIPT:

try: arcpy
except: import arcpy
import itertools
from itertools import combinations

mxd = arcpy.mapping.MapDocument("CURRENT")
layer = arcpy.mapping.ListLayers(mxd)[0]
t0=time.time()

#lines to change
meanV=849
stDev=249
n=3
highLimit=(meanV+stDev/2)*n
arcpy.AddMessage('High limit=%s' %highLimit)

table=[]
for row in arcpy.da.SearchCursor(layer,("Shape@","GRID_CODE")):
    table.append(row)
N=len(table)
count=math.factorial(N)/(math.factorial(n)*math.factorial(N-n))


arcpy.SetProgressor("step", "", 0, count,1)    
dMax=0
for fivePoints in combinations(range(N),n):
    arcpy.SetProgressorPosition()    
    sCur=0
    for i in range(n):
        v=table[fivePoints[i]][1]
        sCur+=v
    if sCur<highLimit:continue
    dCur=0
    for pair in combinations(range(n),2):
        i=fivePoints[pair[0]];p1=table[i][0]
        i=fivePoints[pair[1]];p2=table[i][0]
        dCur+=p1.distanceTo(p2)
    if dCur<dMax:continue
    dMax=dCur;bestPoints=fivePoints
arcpy.AddMessage(bestPoints)
arcpy.AddMessage('Write them done and use "FID" in (...) to select')
arcpy.SelectLayerByAttribute_management(layer, "NEW_SELECTION", '"FID" in %s' %str(bestPoints))

This is very interesting task, perhaps the reason behind down votes is lack of efforts shown by you. Thus try to understand script and make it work.

UPDATE:

• This is how to create toolbox
• This is how to add script to it. Note this script tool has no parameters, just do what I said re location of points in table of content
• and, finally, you wasted your time installing Python. ArcGIS comes with it.

I leave it to you to figure out how to run this tool, after all

Answer 2

I don't feel comfortable answering until you have at least provided some code, but here are some related questions from overflow and here that may get you started.

You need to specify what program or programming language(s) you are comfortable with, so more answers will be relevant.

In ArcGIS the first step would be to convert your rasters to points and use the "Point Distance Tool". You could then find the max distance for each point and iterate through them to maximize the distance between all five points.

Python: Performing a 'far analysis' using ArcGIS for Desktop?

https://stackoverflow.com/questions/24473503/what-is-the-most-efficient-way-to-calculate-the-maximum-distance-of-two-points-i

https://stackoverflow.com/questions/31667070/max-distance-between-2-points-in-a-data-set-and-identifying-the-points