Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a pretty big point feature class in a file geodatabase (~4 000 000 records). This is a regular grid of points with a 100m resolution.

I need to perform a kind of generalization on this layer. For this, I create a new grid where each point lies in the middle of 4 "old" points:

 *     *     *     *
    o     o     o
 *     *     *     *
    o     o     o
 *     *     *     *

[*] = point of the original grid - [o] = point of the new grid

The attribute value of each new point is calculated based on the weighted values of its 4 neighbors in the old grid. I thus loop on all the points of my new grid and, for each of them, I loop on all the points of my old grid, in order to find the neighbors (by comparing the values of X and Y in the attribute table). Once 4 neighbors have been found, we get out of the loop.

There is no methodological complexity here but my problem is that, based on my first tests, this script will last for weeks to complete...

Do you see any possibility to make it more efficient? A few ideas on the top of my head:

  • Index the fields X and Y => I did that but didn't notice any significant performance change
  • Do a spatial query to find the neighbors rather than an attribute-based one. Would that actually help? What spatial function in ArcGIS should do the job? I doubt that, e.g., buffering each new point will prove more efficient
  • Transform the feature class into a NumPy Array. Would that help? I haven't worked a lot with NumPy so far and I wouldn't like to dive into it unless someone tells me it might really help reducing the processing time
  • Anything else?

Thanks in advance!

share|improve this question
    
What version of Arcmap are you using? –  Martin Jan 17 at 13:52
    
Have you considered PostGIS? Is that an option? –  Chad Cooper Jan 17 at 13:56
    
Sorry that I forgot that: ArcGIS 10.1 // Python 2.7 –  Stéphane Henriod Jan 17 at 13:56
    
Nope, PostGIS is unfortunately not an option, my hands are unfortunately quite tied here... At best I can use Oracle with the SDE functions –  Stéphane Henriod Jan 17 at 13:57
add comment

5 Answers 5

For a regular grid, it should be by far more efficient to work in a raster format. Convert your first grid into a raster, the you can resample at the same resolution using a bilinear interpolator but shifting your output image by 1/2 pixel in X and Y, and back again to points if you still need to have points.

EDIT : for complex decisions rules, you can convert each of the fields that you need as a new raster band, then you make four copies of those bands and you shift you raster in the 4 directions by 1/2 pixel (+50, -50), (+50,+50), (-50,-50) and (-50,+50). Then you can use regular map algebra

share|improve this answer
    
Thanks I have actually thought of this solution but I am not sure if / how I can implement the calculation of the new value if in raster format. Let me explain: for each new point (or new raster cell) I need calculate its value as such: I take the value of each of its neighbors. Each of those values has a probability of giving a particular value to the new point. E.g., if one neighbor has the value 202, then it will give the value 3 (with a weight of 1) or the value 11 (with a weight of 5). We then sum up for all the 4 neighbors and find the new value... Not sure whether this is very clear... –  Stéphane Henriod Jan 17 at 14:12
    
PS: the calculation to find the new value can, in some cases, be based on 2 attributes, not just one, which might discard the Raster approach –  Stéphane Henriod Jan 17 at 14:17
    
for your weighted sum, you just need two rasters : one where you resample the product of the weights and the values, the second where you resample the weights only. If you divide the first by the second, you obtain your weighted sum. –  radouxju Jan 17 at 19:46
1  
@StéphaneHenriod - as a suggestion, you might consider editing the question to add these additional specifications. Given the initial question, I think this answer makes tons of sense, but with this new information, Barbarossa's answer looks good. –  nicksan Jan 22 at 1:23
add comment

What if you fed the points into a numpy array and used a scipy cKDTree to look for neighbors. I process LiDAR point clouds with large numbers of points (> 20 million) in several MINUTES using this technique. There is documentation here for kdtree and here for numpy conversion. Basically, you read the x,y into an array, and iterate over each point in the array finding indices of points within a certain distance (neighborhood) of each point. You can use these indices to then calculate other attributes.

share|improve this answer
    
this answer is better than mine –  radouxju Jan 17 at 15:12
    
I like this idea but I don't have scipy on the workstation I'm working on (and no admin rights). If I manage to get this package installed, then I'll give it a try –  Stéphane Henriod Jan 17 at 16:28
add comment

I am with Barbarossa... arcpy cursors are insanely lame, so I only use them to traverse a table or feature class exactly one time. If I can't get the job done in one cycle, I use the cursor to fill up some other kind of data structure and work with that.

If you do not want to hassle with numpy just make a simple python dictionary where you use your coordinates as a simple text key, and fill in the attributes you need for calculation into a list as the value of the dictionary item.

In a second step you can easily get the values you need to calculate a point by simply getting them from your dictionary (which is incredibly fast, because of the dictionaries hashindex of items).

share|improve this answer
    
I actually like your idea with dictionaries and I just implemented it. It indeed works much better... until I actually write the results with a row.insertRow()... Any idea how I can improve this part as well? –  Stéphane Henriod Jan 17 at 16:28
    
I had an similar problem where i had to select some 10.000 points outof 14 Mio. and then delete it. arcpy.cursors where only able to delete about 1 or 2 points per second (!). so i installed pyodbc module to delete them with a single SQL DELETE Statement in only one second. UPDTATEing over SQL will bring you much improvement, as long as you only want to modify attributes...nevertheless you will have to install additional python modules...but it's worth it. –  Jürgen Zornig Jan 17 at 23:19
add comment

Problem with cursors is that you can cycle through them in one way only and you cant go back. Although not recommended, you can populate the feautres into a structure if you are planning to revisiting them.

If you were able to process your features in a single loop, I suggest enabling recycling. It is a parameter on your search featureclass function that lets python to reuse the memory allocated by old features and make traversing the features in a cursor much faster. You can process your grid 80% faster.

Problem is that you can't enable recycling if you are planning to store retrieved features from a cursor.

share|improve this answer
    
I want to explore this "recycle cursor" subject but cannot find any documentation on the ESRI Help. Do you have a link? The Search Cursor has no a recycle parameter. Select_by_Attribute has no such parameter. I see nothing in ENV. –  klewis Jan 17 at 22:24
    
I wrote an article a while back husseinnasser.com/2009/08/when-to-use-recycling-cursor.html?m=1 –  hnasr Jan 18 at 4:42
    
I don't think "reusing cursors" is available through ArcPy, only with the core Arcobjects. –  klewis Jan 20 at 22:25
add comment
up vote 2 down vote accepted

Thanks everybody for your help!

I finally found a very non-pythonic way to solve this issue... What was actually taking the most computing time was to find the 4 neighbors of each point. Rather than using the X and Y attributes (either with an arcpy cursor or within another data structure, such as a python ditionary), I ended up using the ArcGIS tool Generate near table. I assume this takes advantage of the spatial indexes and the performances are obviously much much higher, without me having to implement the index myself.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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