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

Nikos Alexandris nikos.alexandris at felis.uni-freiburg.de
Thu Aug 12 00:56:04 EDT 2010


So here you go using data from spearfish60 [1][2][3]

# in grass64 using v.what.vect (makes no difference, right?)
time v.what.vect --v vector=pareto_ref_points___pareto_ref_100m 
qvect=pareto_grid___OVER___ref_coi_points___pareto_grid___OVER___ref_100m___ON___classification_V___classification_1_200m 
column=gridcell_1 qcolumn=cat

Finding nearest lines...
Finding nearest areas...
26335 categories read from the map
26335 categories exist in the table
26335 categories read from the map exist in the table
26335 records updated
v.distance complete.

real    3m3.687s
user    2m38.910s
sys     0m2.400s

# in grass70 before the latest update
time v.distance --v from=pareto_ref_points___pareto_ref_100m 
to=pareto_grid___OVER___ref_coi_points___pareto_grid___OVER___ref_100m___ON___classification_V___classification_200m_1 
column=gridcell_1 to_column=cat upload=to_attr dmax=0.0

Finding nearest feature...
Finding nearest areas...
26335 categories read from the map
26335 categories exist in the table
26335 categories read from the map exist in the table
26335 records updated
v.distance complete.

real    3m4.270s
user    2m43.140s
sys     0m2.870s

# in grass70 after the latest update ( I only did "make" within the 
/vector/v.distance directory - correct?)

Finding nearest feature...
Finding nearest areas...
Update database...
26335 categories read from the map
26335 categories exist in the table
26335 categories read from the map exist in the table
26335 records updated
v.distance complete.

real    3m5.339s
user    2m42.780s
sys     0m3.040s


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

:-?

Nikos
---

[1] map #1
v.info -t pareto_ref_points___pareto_ref_100m

nodes=26335
points=26335
lines=0
boundaries=0
centroids=0
areas=0
islands=0
primitives=26335
map3d=0

[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

[3] The above maps where created by using data from spearfish60 as follows:

# enter spearfish60
grass64 /spearfish60/user1/

#
g.copy rast=landcover.30m,landcover.30m
g.region rast=landcover.30m -pa


# suppose the landcover.30m is the "source" map
# derive a "reference" high resolution dichotomic map (=100m)
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 radius=1.1 metric=manhattan


# create the 100m maps (source + reference derived from the source)
g.region rast=landcover.30m res=100 -pa
r.mapcalc "landcover.100m = landcover.30m"
r.mapcalc "pareto_ref_100m = pareto_ref"
g.remove rast=pareto_ref


# create dichotomic rangeland maps supposedly being classification results
g.region rast=landcover.30m -pa
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())"


# suppose those maps come from data of low resolution (finally =200m)
g.region rast=landcover.30m res=200 -pa
r.grow input=coi_1 output=pareto_classification_200m_1 radius=1.1 
metric=manhattan
r.grow input=coi_2 output=pareto_classification_200m_2 radius=1.1 
metric=manhattan
g.mremove rast=coi_[12] -f

# below using the scripts that can be found at ticket 
# pareto step 1
pareto_1_vectorise_rasters.py reference_raster=landcover.100m 
reference_coi_rasters=pareto_ref_100m 
classification_rasters=pareto_classification_200m_1,pareto_classification_200m_2

[...]
   Vectorised classifications maps for next step:
[  'pareto_ref_points___pareto_ref_100m',
   'pareto_ref_coi_points___pareto_ref_100m',
   'pareto_classification_V___pareto_classification_1_200m',
   'pareto_classification_V___pareto_classification_2_200m']

# pareto step 2 (this also takes time!)
pareto_2_create_lowres_vector_grid.py highres=100 lowres=200 --v


# pareto step 3 (timing v.what.vect here!)
time v.what.vect --v vector=pareto_ref_points___pareto_ref_100m 
qvect=pareto_grid___OVER___ref_coi_points___pareto_grid___OVER___ref_100m___ON___classification_V___classification_1_200m 
column=gridcell_1 qcolumn=cat

All above steps repeated within grass70 before the update (should make no 
difference, right?) and only the last script execution (pareto_3) skipped, 
instead the v.distance command tested directly.

!!!
a. my script does not work in grass70
b. v.what.vect seems to be broken after latest update?

--> just test "v.distance"
!!!

time v.distance --v from=pareto_ref_points___pareto_ref_100m 
to=pareto_grid___OVER___ref_coi_points___pareto_grid___OVER___ref_100m___ON___classification_V___classification_200m_1 
column=gridcell_1 to_column=cat upload=to_attr dmax=0.0


More information about the grass-dev mailing list