[GRASS-user] how to measure distances in a constrained space

GRASS GRASS user
1

Hello,

I would like to measure distances between points in a constrained space, such as distances between cities, but only taking into account distances by land, and so excluding any water bodies. Concretely, this means finding those cities in a first set of cities which are at least 150km from the nearest city in the a second set of cities.

At this stage, I only see the possibility of creating a cost surface with cell values on land masses equal to the resolution and null cells on water masses and then using r.cost with null_cost set to a very high level and start_points to the second set of cities, then v.what.rast to sample the result at the location of the cities of the first set. Does that sound reasonable ? At the moment I'm getting very weird values, but I'll have to double-check my operations.

Does anyone see a way of doing this with vector operations ?

Moritz

2

I've implemented exactly the same idea as yours in this script:
http://www.geeitema.org/doc/guenmap//docs/v.costdist.mat.zip

...with the particularity that it stores the computed distances in the attribute table of the second set of cities, which was needed for my application.
You can try it, and see if it works for you.

Regards!
ƒacu.-

#
# Description:
# v.costdist.mat determines de cumulative cost of moving between each pair of points
# from two given maps over a cost surface
# e.g.
# - obtaining the "true" distance matrix of a set of noise measures in a city
# considering buildings as non-transparent areas for sound
#
# Notes:
# Results are stored as new columns in "to" map's attributes tables
# Try and keep the number of input sites to under a few dozen, as the
# process is very computationally expensive. You might consider creating
# a few hundred MB RAM-disk for a temporary mapset to avoid thrashing your
# hard drive (r.cost is heavy on disk IO).

Moritz Lennert escribió:

Hello,

I would like to measure distances between points in a constrained space, such as distances between cities, but only taking into account distances by land, and so excluding any water bodies. Concretely, this means finding those cities in a first set of cities which are at least 150km from the nearest city in the a second set of cities.

At this stage, I only see the possibility of creating a cost surface with cell values on land masses equal to the resolution and null cells on water masses and then using r.cost with null_cost set to a very high level and start_points to the second set of cities, then v.what.rast to sample the result at the location of the cities of the first set. Does that sound reasonable ? At the moment I'm getting very weird values, but I'll have to double-check my operations.

Does anyone see a way of doing this with vector operations ?

Moritz

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

3

Moritz Lennert escribió:

Hello,

I would like to measure distances between points in a constrained space, such as distances between cities, but only taking into account distances by land, and so excluding any water bodies. Concretely, this means finding those cities in a first set of cities which are at least 150km from the nearest city in the a second set of cities.

On 12/11/08 19:46, Facundo Muñoz wrote:
> I've implemented exactly the same idea as yours in this script:
> http://www.geeitema.org/doc/guenmap//docs/v.costdist.mat.zip
>
> ...with the particularity that it stores the computed distances in the
> attribute table of the second set of cities, which was needed for my
> application.
> You can try it, and see if it works for you.

Thank you ! It is actually a bit overkill for me as I don't need the distance from every city in set 1 to every city in set 2, but only the distance for every city in set 1 to the closest city in set 2. As I have over 1000 cities in one and 556 cities in the other set, it is a bit heavy for your script.

I actually now get a correct result with the following steps (my original cost map was wrong, that's why my results where weird):

- create a cost map where all land parts are equal to the resolution and all water masses are null
- run r.cost with start_points=set2 (actually no need for null_cost) tio create the distance map
- run v.what.rast vector=set1 raster=distance column=dist_column

And results seem perfectly correct.

Now just for curiosity, any possible solution in vector space ? I thought about v.net.visibility, but it would need to be altered to allow a visibility graph from vector1 to vector2, instead of only within vector1. Then I guess one could use something like v.net.iso or loop through v.net.path...

Moritz

4

Moritz Lennert wrote:

>> I would like to measure distances between points
>> in a constrained space, such as distances between cities,
>> but only taking into account distances by land, and so
>> excluding any water bodies. Concretely, this means finding
>> those cities in a first set of cities which are at least
>> 150km from the nearest city in the a second set of cities.

Facundo Muñoz wrote:

> I've implemented exactly the same idea as yours in this script:
http://www.geeitema.org/doc/guenmap//docs/v.costdist.mat.zip
>
> ...with the particularity that it stores the computed distances in the
> attribute table of the second set of cities, which was needed for my
> application.
> You can try it, and see if it works for you.

I actually now get a correct result with the following
steps (my original cost map was wrong, that's why my
results where weird):

- create a cost map where all land parts are equal to the
resolution and all water masses are null
- run r.cost with start_points=set2 (actually no need for
null_cost) tio create the distance map
- run v.what.rast vector=set1 raster=distance column=dist_column

And results seem perfectly correct.

Now just for curiosity, any possible solution in vector
space ? I thought about v.net.visibility, but it would need
to be altered to allow a visibility graph from vector1 to
vector2, instead of only within vector1. Then I guess one
could use something like v.net.iso or loop through
v.net.path...

just some extra notes:
- r.cost in GRASS 6.4 should be 60x faster than in older versions
- have a look at the new r.grow.distance in GRASS 7, even faster
- FWIW when I've done this I've made the cost map 1,NULL and multiplied
  the cost distance by map res. knights move is essential.
- See the v.surf.icw script in wiki addons. it is like IDW interpolation
  but replaces euclidean distance with true (cost) distance and goes
  around corners/coastlines.
- Somewhere I had a script for generating distance matrix tables for
  distance between sites taking obstacles into account using r.cost.
  I'd have to hunt for a backup copy.

It would definitely be interesting to hear if you figure something out
with v.net.visibility. you are looking for as-the-bird-flies distance,
not road distance, right?

see ongoing geodesic distance thread on the PROJ4 list for more on the
dangers of measuring long distances on a grid map projection.

Hamish

5

On 13/11/08 04:14, Hamish wrote:

Moritz Lennert wrote:

I would like to measure distances between points
in a constrained space, such as distances between cities,
but only taking into account distances by land, and so
excluding any water bodies. Concretely, this means finding
those cities in a first set of cities which are at least
150km from the nearest city in the a second set of cities.

Facundo Muñoz wrote:

I've implemented exactly the same idea as yours in this script:

http://www.geeitema.org/doc/guenmap//docs/v.costdist.mat.zip

...with the particularity that it stores the computed distances in the
attribute table of the second set of cities, which was needed for my
application.
You can try it, and see if it works for you.

I actually now get a correct result with the following
steps (my original cost map was wrong, that's why my
results where weird):

- create a cost map where all land parts are equal to the
resolution and all water masses are null
- run r.cost with start_points=set2 (actually no need for
null_cost) tio create the distance map
- run v.what.rast vector=set1 raster=distance column=dist_column

And results seem perfectly correct.

Now just for curiosity, any possible solution in vector
space ? I thought about v.net.visibility, but it would need
to be altered to allow a visibility graph from vector1 to
vector2, instead of only within vector1. Then I guess one
could use something like v.net.iso or loop through
v.net.path...

just some extra notes:
- r.cost in GRASS 6.4 should be 60x faster than in older versions

Am using grass7, can confirm that it is really fast.

- have a look at the new r.grow.distance in GRASS 7, even faster

r.grow.distance creates a distance map from each point, but there is no way of telling it anything about null cells or barriers. Even when setting a mask, r.grow.distance fills in the whole region, which seems like a bug to me.

- FWIW when I've done this I've made the cost map 1,NULL and multiplied
  the cost distance by map res.

What's the difference to directly using a cost map res,NULL ?

knights move is essential.

I guess that's an issue of precision, or ? For me "normal" moves correctly identifies the cities I'm looking for.

- See the v.surf.icw script in wiki addons. it is like IDW interpolation
  but replaces euclidean distance with true (cost) distance and goes
  around corners/coastlines.

Interesting, but:
"Input data points [...] but should be kept within a few dozen as the module becomes very computationally expensive"

IISC, Facundo's module is based on that, or ?

So also not really an option with >1000 points.

- Somewhere I had a script for generating distance matrix tables for
  distance between sites taking obstacles into account using r.cost.
  I'd have to hunt for a backup copy.

Well the above works, i.e. just setting obstacles to NULL in the cost map. Obviously this only works if obstacles represent total barriers, not just increased friction.

It would definitely be interesting to hear if you figure something out
with v.net.visibility. you are looking for as-the-bird-flies distance,
not road distance, right?

Road distance would be even better, but the v.net.* modules unfortunately don't scale at all...

see ongoing geodesic distance thread on the PROJ4 list for more on the
dangers of measuring long distances on a grid map projection.

I'll have a look, thanks !

Moritz

6

Moritz Lennert escribió:

- See the v.surf.icw script in wiki addons. it is like IDW interpolation
but replaces euclidean distance with true (cost) distance and goes
around corners/coastlines.

Interesting, but:
"Input data points [...] but should be kept within a few dozen as the module becomes very computationally expensive"

IISC, Facundo's module is based on that, or ?

Ja! you got me!
No, the module is not based on v.surf.icw. I just used it as a "template" or a guide for programming mine, since i had never wrote one. And since I also make intensive use of r.cost, I just kept that phrase in the description.

Regards!
ƒacu.-

7

Moritz Lennert wrote:

> - have a look at the new r.grow.distance in GRASS 7, even faster

r.grow.distance creates a distance map from each point, but there is no
way of telling it anything about null cells or barriers. Even when
setting a mask, r.grow.distance fills in the whole region, which seems
like a bug to me.

A mask just forces some input cells to be null. As a general
principle, modules don't distinguish between cells which are null in
the map itself and those which are null due to the mask.

In any case, the specific purpose of r.grow.distance is:

   r.grow.distance generates a raster map representing the distance to
   the nearest non-null cell in the input map.

This is its only purpose, and it tries to do it as efficiently as
possible. In particular, it operates row-by-row, so it doesn't require
the map to fit into memory, nor does it use a segment cache (it does
use a temporary file, but this is read linearly).

There's no way to make the algorithm handle barriers, as it relies
upon the fact that the distance from a point increases monotonically
with each row.

Also, r.grow.distance will give true circles (for metric=euclidean and
metric=squared) rather than a polygonal approximation.

However, I have just noticed a bug: it performs the square root
calculation for all metrics, not just euclidean. I'll fix this
shortly.

--
Glynn Clements <glynn@gclements.plus.com>

8

Hamish escribió:

- See the v.surf.icw script in wiki addons. it is like IDW interpolation
but replaces euclidean distance with true (cost) distance and goes
around corners/coastlines.

Moritz:

Interesting, but:
"Input data points [...] but should be kept within a few dozen as the
module becomes very computationally expensive"

Partial weights are constructed by running r.cost starting from each
input point. A much faster r.cost will mean it is practical to run it
with many more points, but disk space/IO may still be a problem.
(mapset on ramdisk solution)

I will have some new data to process with it in the next weeks, but
perhaps ~200 points might be ok now. Of course resolution matters as much
as number of starting points, I speak of 2000x2000 or so raster region
sizes.

Hamish

9

Moritz:

Well the above works, i.e. just setting obstacles to NULL
in the cost map. Obviously this only works if obstacles
represent total barriers, not just increased friction.

if cost map is all "1" and you multiply by map resolution later to get
distance, then you could simulate increased friction by having some cost
cells as 2,3,4, etc... twice as slow is effectively twice as far.

I am working in waterways, I have wondered about generating an upstream/
downstream cost map to use with v.surf.icw so interpolation following
river bends prefers not to mix pollutants/invasive species from tributaries
upstream as readily as downstream... ie make it more expensive to move
upstream, but not impossible (tidal exchange, diffusion, turbulence, etc).
That module does not care about absolute distance, just that all use the
same method.

Hamish