Hello everyone,
Time for a little feedback I think. First of all, thanks to both of you for the incredible help.
And more, thank you for the invaluable good practice tips, I had no idea of where to find such advice, still not but I'm keeping all you said in mind.
As for answering your questions quickly:
1) I am using arcgis 10.2.2 and had no idea there was a new cursor: I'll look at the correct version web help now.
2) The variable 'contour' is the shapefile "perimeter.shp" descibing the study area form, which I forgot to rename in the process of trying to make the code a bit clearer for posting.
5) "VALUE > 180002": this was the value at which the precedent run had crashed, so I extracted the points that remained to be processed before launching a new run. Dirty, but in 5-10 times I had them all. Still, I hope I'll learn enough to avoid that problem another time.
3), 4) and 6): thanks for the advice
As I was really running out of time and courage, I have implemented the last part of the code in ModelBuilder, and I don't think I'll be able of doing it in Arcpy right now. So I'm posting the result in order for others to look at if they need ideas.
It aims at extracting the pixel distribution histogram for a great number of zones. In the picture below you can see how it works:
In the background the land cover map, resolution 15m x 15m: I wanted to extract the land cover in a 150m radius around each pixel.
For this I computed discs covering that area, and to go faster I made fishnets made of disks instead of single discs
As each disks covers an area of 21 x 21 = 441 pixels, and Zonal Histogram won't handle overlapping features, I had to make as many disk-fishnets to be able to cover the whole number of pixels. (two fishnets are shown in the picture)
All this runs in the arcpy code provided some patience and re-starting the code after crashes.
Then I could extract the land cover distribution for each disc-fishnet, in model builder using a double loop (i.e. two nested models).
The 441 fishnet files are stored in the file zones.gdb. They have 3 relevant fields:
1) "orig_id" gives a unique value for each disk in order to delimit the zones for the "Zonal histogram" tool.
2) "Id2" = Int(FID/1000) +1 allows to group disks in 46 lots of 1,000.
3) A third field of 0 and 1 allowed to filter out disks which have their centers outside the land cover raster
The model has basically 3 steps:
4) Main model: for each feature class (fishnet), I convert the unique-value-per-disk field, "orig_id", to raster.
5) Nested model: another loop extracts a fragment of this raster by groups of 1,000 using the ID2 field. For each value of Id2, I make a feature layer of only the features with Id2 = for-loop- increment*, and extract the fragment from the whole fishnet fishnet with "extract by attributes".
6) Nested model: each fragment of orig_id.img is used as a zone field for zonal histogram.
Main loop:
Nested loop:
Here are the main problems I encountered and the work arounds I found:
2) Zonal histogram won't run over more than 1,000 features at a time, so I needed to split the job. the "+1" is because in the nested model, which is a for loop, iteration over the value 0 always fails, (Error message: "the precondition is false" and then in the downloop steps "the inputs are not current") so I looped over 0:47 instead of 0:46 and the loop bug in 0 does not count anymore.
3) if you leave them the Zonal Histogram tool crashes completely.
4) Back when I was trying arcpy and using feature classes as zone input for Zonal Histogram, it would automatically convert them to raster and automatically name them, with no workaround. They piled up in the workspace and made the loop eventually crash. By converting the fragment to raster this is avoided.
5) "Zonal histogram" will ignore an output of "select layer by attributes", so "extract" is compulsory. Also, the SQL query within "make feature layer" does not work so I had to use both this and "select layer by attributes".
All in all, the loop takes 1 hour to loop over 1 fishnet. Limiting steps are 'extract by mask' (15s) and 'zonal histogram' (1mn10s). The expected calculation time for the 441 fishnets will be 18 days, but with parallelisation in 10 simultaneous runs this can be brought down to 2 days.
And this is better than the 2 years with which I started I just hope the loop will be stable over 40 feature classes.