1

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)
...

enter image description here

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`
8
  • 2
    Please edit the question to contain the code you have developed to date. It should include the names(s) of data sources and the version of ArcGIS in use.
    – Vince
    May 16, 2016 at 3:11
  • 1
    A picture of what your raster looks like might help. There is a function called 'Euclidean Distance' but I don't know if it will help in this case.. you might need to extract by attributes the values you want, convert to point and use a tool like Point Distance to find the distances between each cell. May 16, 2016 at 3:19
  • Question edited. I cannot provide the code developed to date as I am stuck in this optimization problem and don't even know how to start it with code, as it's the first time doing so.
    – tom
    May 16, 2016 at 23:13
  • You have provided R code but this is tagged as ArcGIS?
    – GISHuman
    May 18, 2016 at 17:54
  • 1
    Your problem has no solution. Redefine 2nd condition as sum being above certain threshold and it will
    – FelixIP
    May 21, 2016 at 20:54

2 Answers 2

3
+50

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:

enter image description here

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

5
  • Copy it to text file and call it points.py. Create toolbar and right click, script
    – FelixIP
    May 22, 2016 at 23:44
  • Much appreciated. I am following your instructions and copying the code in the Python window of ArcMap 10.2, while changing: mxd = arcpy.mapping.MapDocument("G:\ArcGIS\points_file.mxd") meanV=custom according to the initial raster file stDev=custom according to the initial raster file I get the error (see initial post now edited), which I believe is related to SetProgressor. Should I previously download itertools? Or is it embedded in ArcMap 10.2? I have been searching and by writing "import itertools" should have been enough.
    – tom
    May 22, 2016 at 23:52
  • I have no clue how to work such long script from python window, see help on how attach it to tool
    – FelixIP
    May 22, 2016 at 23:54
  • I have installed PythonWin 2.7.5 win32 and when I run it, it says "Exception raised while running script". I am also trying to run it the Python window in ArcMap 10.2 and gives the error with SetProgressor shown above. I have Python 27 and Python 35, how else should I try to run the script? I've google about creating a toolbar to run the script, but still I've not figured out how. Could you please tell me your versions of Python and how exactly you run it, so I will emulate it? Thanks
    – tom
    May 24, 2016 at 2:42
  • See update. It is a tool to run FROM mxd.
    – FelixIP
    May 24, 2016 at 5:29
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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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