[GRASS-dev] About v.distance, v.what.vect (wrt "count points within...").

[Nikos Alexandris]

Nikos A:

(MM & ML, I really appreciate the time you spend on this)

[...]
 
> AFAIU, MarkusM's change only affects dmax>0, so it's normal that you
> don't see a difference.

Yep, it's clear.

> > [2] map #2 v.info -t pareto_...<Very Long Name>

> > nodes=26956 points=0 lines=0 boundaries=26955 centroids=6650
> > areas=6650 islands=1 primitives=33605 map3d=0
 
> So, you have 6000+ areas times 26000+ points.

Well, this is just a "spearfish" small example. Worst case beforehand here is: 
areas=2376  _vs_  points=593273 (in a region of rows:1201, cols:550, 
cells:660550 -- but this does not make any difference though)

> > Hmmm... ? So, if use my 600000+ points maps, you can imagine why I
> > report timings>20h.

> You mean that it's normal that it takes so long, or that it is not
> normal, and that's why we can imagine that you report it ? ;-)

It was a (very) wild guess, smthn like "3m5.339s for 23.000, how much for 600K 
thinking that time-consumation does _not_ increase linearly, but more 
aggressive).

> Here's an attempt with a 100x100 grid (=10000 areas) and the 600000
> points randomly created:
> 
> v.mkgrid mygrid n=100,100 --o

Can someone confirm that even this takes (a bit) long (or is it only my box)?

> time v.distance from=mypoints at sqlite to=mygrid upload=to_attr
> column=cat_municip to_column=cat dmax=0.0

> real	2m22.505s

GRASS 6.4.0svn (nc_spm_08):~ > real    1m42.860s  [ :-) | :-( ]

(shrug)

Don't know what's wrong here!? To make the long story short following my 
scripts:

1. pareto_1 converts all pixels of Landsat7 samples to (centroids ->) points. 
Yes, res=30m, many pixels (e.g. the 593273 above). Includes only: g.region, 
r.to.vect, v.db.dropcol.

2. pareto_2 loops over the point(s) maps and creates "low-resolution" vector 
grid (using v.mkgrid, e.g.  2376 boxes that correspond to a raster with 
res=500m). Includes only: g.region, v.in.region, v.select, g.remove, v.mkgrid

3. pareto_3 just counts points in polygons. Includes: v.db.addcol, v.what.vect 
(this is the time-killer), v.support, v.db.addcol, v.db.update (also 
inefficiently executed but not the real problem)

(...the rest of pareto_? is another story...)


I will give another go with the spearfish60 dataset but this time keeping the 
high-resolution to 30m and the low resolution to 500m. In detail, all of the 
commands that have to be executed are:

--%<---
grass64 /spearfish60/user1/

g.copy rast=landcover.30m,landcover.30m

g.region rast=landcover.30m -pa # cells: 294978

r.mapcalc "rangeland.30m = if((landcover.30m == 51 || landcover.30m == 71 || 
landcover.30m == 81 || landcover.30m == 92), 2, null())"

r.grow input=rangeland.30m output=pareto_ref_30m radius=1.1 metric=manhattan

r.mapcalc "coi_1 = if((landcover.30m == 51 || landcover.30m == 71 || 
landcover.30m == 81 ), 2, null())"

r.mapcalc "coi_2 = if(( landcover.30m == 71 || landcover.30m == 92), 2, 
null())"

g.region rast=landcover.30m res=500 -pa

r.grow input=coi_1 output=pareto_classification_500m_1 radius=1.1 
metric=manhattan

r.grow input=coi_2 output=pareto_classification_500m_2 radius=1.1 
metric=manhattan

g.mremove rast=coi_[12] -f

pareto_1_vectorise_rasters.py reference_raster=landcover.30m 
reference_coi_rasters=pareto_ref_30m 
classification_rasters=pareto_classification_500m_1,pareto_classification_500m_2

pareto_2_create_lowres_vector_grid.py highres=30 lowres=500 --v


# here would go the "pareto_3" script, running many times v.what.vect 
depending on the number of "area-maps" and naturally the "poin-maps". However, 
one test-command is enough to demonstrate:

time v.what.vect --v vector=pareto_ref_points___pareto_ref_30m 
qvect=pareto_grid___OVER___ref_coi_points___pareto_grid___OVER___ref_30m___ON___classification_V___classification_1_500m 
column=gridcell_1 qcolumn=cat
-->%---

Now this is crawling and I am waiting [*]. Didn't wanna let you waiting too 
long to know that this keeps me waiting - so I've hit the Send button :-p

Nikos
---
[*] You Waited Too Long, Casey Bill Weldon (1936)