[GRASS-user] [Bulk] Re: Two simple questions

Hi Luciano:
(Returning to the maillist - maybe this will help others...)

On 05/11/2012 18:29, Luciano La Sala wrote:

Two simple questions

Thank you Micha, I found the icons right under my nose.

To tell you the truth I am kind of lost here. I am trying to do the interpolation from within QGIS using the raster map of risk I did following your instructions (JPEG attached).

Pixel values, as calculated using the “r.mapcalculator” module in GRASS, are (from the lightest yellow to the darkest red):

0.0406

0.1422

0.4268

0.5284

0.7154

0.8983

I first used the module “r.surf.idw” module (Screen_01) which runs successfully but produces an entirely black raster map as output (Screen_02). All the pixels of this raster contain a value of 1.

There is obviously something wrong here. Any tips? Should I use “v.surf.idw” instead?

I found on an old maillist post that r.surf.idw works only with integer values. That's why you're getting everything "rounded" to 1. To work around this, you can try two options:

* revert back to vectors, as you suggest. (v.surf.idw accepts decimal values)
     r.to.vect -z in=risk_raster out=risk_vector feature=point
     v.surf.idw in=risk_vector out=risk_idw layer=0

OR

* multiply your risk raster by some constant, say 1000, to get integer values. Then run r.surf.idw on that, and divide the result by the same 1000 again to get back to your normalized values.
     r.mapcalc risk_integer=risk_raster*1000
     r.surf.idw in=risk_integer out=risk_idw_tmp
     r.mapcalc risk_idw=risk_idw_tmp/1000.0
     # Use the decimal value 1000.0 to insure that the result is floating point.

Regards,

Micha

Thank you so very much for your time Micha.

Luciano

*De:*Micha Silver [mailto:micha@arava.co.il]
*Enviado el:* Monday, November 05, 2012 4:01 AM
*Para:* Luciano La Sala
*Asunto:* [Bulk] Re: Two simple questions

On 05/11/2012 04:42, Luciano La Sala wrote:

    Hello Micha,

    Where do you access the
    modules///v.surf.idw/and///r.surf.idw/from?I just can’t seem to
    find them!

If you run GRASS on its own (not the QGIS plugin), these modules are under the "Raster->Interpolate Surfaces"menu. In addition, all GRASS commands can be run straight from the command line.

    Thanks!

    Luciano

    This mail was received via Mail-SeCure System.

--
Micha Silver
GIS Consulting
052-3665918
http://www.surfaces.co.il

This mail was received via Mail-SeCure System.

--
Micha Silver
GIS Consulting
052-3665918
http://www.surfaces.co.il

···

Hi Luciano:

Keeping the discussion on the grass-users list

On 11/7/2012 10:57 PM, Luciano La Sala wrote:

I tried both options (your message below); i.e. reverting back to vectors and multiplying my risk raster by a constant (in this case 1000) to get integer

values which “r.surf.idw” is capable of dealing with. Option 1 does not work, and I don’t know why. Option 2 seemed more intuitive and simple…

(clipped)

Then run r.surf.idw on that and I get Screen_03, the weirdest map ever!

I think you’ve done everything right. One detail you might not be aware of: the r.surf.* modules will always work in the full “current computational region”. So interpolation continues, as best as it can, right to the edge of the region, creating those funny looking angles and spikes outside of argentina.

The way to do this is create a mask raster. In GRASS you can limit almost all raster command to any arbitrary area by using a mask. So:

1- take a vector polygon of the whole country, and import it into GRASS.

2- Then use “v.to.rast argentina_poly type=area use=val value=1 out=argentina_rast” to convert it to a raster.

3- Next set this argentina raster as the MASK with “r.mask argentina_rast”

4- and now rerun the interpolation

Finally, I divide the result by 1000.0 again (not 1000, like you said, to insure that the result is a floating point) to get back to my normalized values, but the same weird raster map remains (Screen_04).

Check the values you got with r.info.

Where is it that I am going so wrong? To do all this I didn’t use QGIS, so I am not referring the message to the list.

Very best,

Luciano

<<…>> <<…>> <<…>> <<…>>

-----Mensaje original-----
De: Micha Silver [mailto:micha@arava.co.il]
Enviado el: Tuesday, November 06, 2012 6:20 AM
Para: Luciano La Sala
CC: grass-user
Asunto: [Bulk] Re: [Bulk] Re: Two simple questions

Hi Luciano:

(Returning to the maillist - maybe this will help others…)

On 05/11/2012 18:29, Luciano La Sala wrote:

Two simple questions

Thank you Micha, I found the icons right under my nose.

To tell you the truth I am kind of lost here. I am trying to do the

interpolation from within QGIS using the raster map of risk I did

following your instructions (JPEG attached).

Pixel values, as calculated using the “r.mapcalculator” module in

GRASS, are (from the lightest yellow to the darkest red):

0.0406

0.1422

0.4268

0.5284

0.7154

0.8983

I first used the module “r.surf.idw” module (Screen_01) which runs

successfully but produces an entirely black raster map as output

(Screen_02). All the pixels of this raster contain a value of 1.

There is obviously something wrong here. Any tips? Should I use

“v.surf.idw” instead?

I found on an old maillist post that r.surf.idw works only with integer values. That’s why you’re getting everything “rounded” to 1. To work around this, you can try two options:

  • revert back to vectors, as you suggest. (v.surf.idw accepts decimal

values)

r.to.vect -z in=risk_raster out=risk_vector feature=point

v.surf.idw in=risk_vector out=risk_idw layer=0

OR

  • multiply your risk raster by some constant, say 1000, to get integer

values. Then run r.surf.idw on that, and divide the result by the same

1000 again to get back to your normalized values.

r.mapcalc risk_integer=risk_raster*1000

r.surf.idw in=risk_integer out=risk_idw_tmp

r.mapcalc risk_idw=risk_idw_tmp/1000.0

Use the decimal value 1000.0 to insure that the result is

floating point.

Regards,

Micha

Thank you so very much for your time Micha.

Luciano

*De:*Micha Silver [mailto:micha@arava.co.il mailto:micha@arava.co.il]

Enviado el: Monday, November 05, 2012 4:01 AM

Para: Luciano La Sala

Asunto: [Bulk] Re: Two simple questions

On 05/11/2012 04:42, Luciano La Sala wrote:

Hello Micha,

Where do you access the

modules///v.surf.idw/and///r.surf.idw/from?I just can’t seem to

find them!

If you run GRASS on its own (not the QGIS plugin), these modules are

under the "Raster->Interpolate Surfaces"menu. In addition, all GRASS

commands can be run straight from the command line.

Thanks!

Luciano

This mail was received via Mail-SeCure System.

Micha Silver

GIS Consulting

052-3665918

http://www.surfaces.co.il

This mail was received via Mail-SeCure System.

Micha Silver

GIS Consulting

052-3665918

http://www.surfaces.co.il

This mail was received via Mail-SeCure System.

Maybe an example will make things clearer:

https://picasaweb.google.com/115558684812103743786/Maps?authuser=0&feat=directlink

Here are three images that I made as follows:
First I imported point locations (cities and towns, in this case) into three separate GRASS point vectors. THen I buffered each with a different distance. Next I converted the area vector buffers to rasters using a different value (weighting) for each. (I also had a MASK defined to stop interpolation at the country border) I did:
  v.to.rast villages_buff out=villages type=area use=val value=1
  v.to.rast towns_buff out=towns type=area use=val value=3
  v.to.rast cities_buff out=cities type=area use=val value=7

Before running the r.mapcalc , I converted NULL values to zero in each of the maps (otherwise, as we discussed, r.mapcalc returns NULL when any one of the maps is NULL). So

r.null cities null=0
r.null villages null=0
r.null towns null=0

Now a simple :
r.mapcalc risk="cities+villages+towns)/11.0

gives me a weighted risk map (smilar to the one you already have). Have a look at the initial image of the three I posted to Picassa web. Grey areas are no value (NULL). The colored areas have *discrete* values (the legend with continuous coloring is misleading): 100,300,700,800, 1000,1100 based on the weighting I gave, and combinations of those weights where there are overlaps.

Next I did the interpolation as follows :
r.mapcalc risk_1000=risk*1000.0
r.surf.idw --o in=risk_1000 out=risk_idw npoints=500

I chose a fairly large "npoints" value to get some smoothing. The result is the IDW risk map. Note two things-
* All the original values are preserved. That's how IDW works.
* Everywhere that there was NULL in the original weighted map got some interpolated value. In those areas there is a continuous gradient of colors, since each pixel got some average of the 500 nearest points from the interpolation.

Next, just as another example, I ran r.neighbors on the weighted risk map:
r.neighbors --o in=risk_1000 out=risk_neigh size=15

The "size=15" parameter sets a window of 15 X 15 pixels. So each value in the "nearest neighbors" raster gets the average of a 15X15 window from the original. This is NOT an interpolation, just a smoothing filter. This resulted in the third image. Note here that most of the areas that were NULL in the original stayed NULL. And areas that had some value previously got "averaged out".

Maybe this will help you decide what interpolation to use...

On 11/10/2012 01:27 AM, Luciano La Sala wrote:

Micha,

Looking very closely to pixels of both rasters, they are identical meaning that “r.surf.idw” performed on that raster did not produced any smoothing. I also tried playing with the only algorithm parameter available to be changed (number of interpolation points) and setting it to a range of values between 4 and 20. Nothing changed. All the maps look exactly the same.

Do you think this problem is somehow related to the scores by which I multiplied each variable (border, poultry farm, wetland, etc.) in the first place?

(Argentina*1.0 + Airports*2.0 + Border*3.0 + Poultry*4.0 + Wetlands*5.0) / 15

Best luck and thank you for your valuable time and patience.

Luciano

*De:*Micha Silver [mailto:micha@arava.co.il]
*Enviado el:* Thursday, November 08, 2012 3:57 AM
*Para:* Luciano La Sala
*Asunto:* [Bulk] Re: Still in trouble

On 08/11/2012 01:45, Luciano La Sala wrote:

    Thank you Micha,

    I will subscribe to the GRASS-users list later tonight.

    You raised an interesting point about using a mask. I will keep
    that in mind. However, if you take a look at Screen_2 (raster
    multiplied by 1000) and Screen_3 (“r.surf.idw” performed on that
    raster), they look exactly the same except for the funny area
    outside of Argentina’s border.

Read up on the Inverse Distance Weighted algorithm. It gives a very high weight to pixels that are very close, and a much smaller weight as the distance from the pixels gets larger. So I wouldn't expect much of a change using that interpolation algorithm. If you look closely at some of the areas in the north center of the map I think you can see areas where the points got "smoothed"out. Yuo can play with the algorithm parameters to get more smoothing of the results.

    Having said that (correct me if I am wrong here), I believe that
    the interpolation function did not work. If it had, shouldn’t we
    be looking at a new raster map (on Screen_3) which is different
    from the one to which we applied the “r.surf.idw” module (Screen_2)?

    Any ideas? Is is enough to multiply my original raster by 1000? As
    I showed in my previous email, that doesn’t yield to integer values.

    Thanks again! Lucho

    *De:*Micha Silver [mailto:micha@arava.co.il]
    *Enviado el:* Wednesday, November 07, 2012 7:19 PM
    *Para:* Luciano La Sala
    *CC:* GRASS user List
    *Asunto:* Re: Still in trouble

    Hi Luciano:

    Keeping the discussion on the grass-users list

    On 11/7/2012 10:57 PM, Luciano La Sala wrote:

        I tried both options (your message below); i.e. reverting back
        to vectors andmultiplyingmyrisk raster byaconstant (in this
        case1000) to get integer

        values which "r.surf.idw" is capable of dealing with. Option 1
        does not work, and I don't know why. Option 2 seemed more
        intuitive and simple...

    (clipped)

    Then run r.surf.idw on that andI get Screen_03, the weirdest map ever!

    I think you've done everything right. One detail you might not be
    aware of: the r.surf.* modules will always work in the full
    "current computational region". So interpolation continues, as
    best as it can, right to the edge of the region, creating those
    funny looking angles and spikes outside of argentina.

    The way to do this is create a mask raster. In GRASS you can limit
    almost all raster command to any arbitrary area by using a mask. So:

    1- take a vector polygon of the whole country, and import it into
    GRASS.

    2- Then use "v.to.rast argentina_poly type=area use=val value=1
    out=argentina_rast" to convert it to a raster.

    3- Next set this argentina raster as the MASK with "r.mask
    argentina_rast"

    4- and now rerun the interpolation

        Finally, Idivide the result by 1000.0 again (not 1000, like
        you said,to insure that the result isafloating point)to get
        back to my normalized values, but the same weird raster map
        remains (Screen_04).

    Check the values you got with r.info.

        Whereis it that I am going so wrong? To do all this I didn’t
        use QGIS, so I am not referring the message tothe list.

        Very best,

        Luciano

        <<...>> <<...>> <<...>> <<...>>

        -----Mensaje original-----
        De: Micha Silver [mailto:micha@arava.co.il]
        Enviado el: Tuesday, November 06, 2012 6:20 AM
        Para: Luciano La Sala
        CC: grass-user
        Asunto: [Bulk] Re: [Bulk] Re: Two simple questions

        Hi Luciano:

        (Returning to the maillist - maybe this will help others...)

        On 05/11/2012 18:29, Luciano La Sala wrote:

        > Two simple questions

        >

        > Thank you Micha, I found the icons right under my nose.

        >

        > To tell you the truth I am kind of lost here. I am trying to
        do the

        > interpolation from within QGIS using the raster map of risk I
        did

        > following your instructions (JPEG attached).

        >

        > Pixel values, as calculated using the “r.mapcalculator”
        module in

        > GRASS, are (from the lightest yellow to the darkest red):

        >

        > 0.0406

        >

        > 0.1422

        >

        > 0.4268

        >

        > 0.5284

        >

        > 0.7154

        >

        > 0.8983

        >

        > I first used the module “r.surf.idw” module (Screen_01) which
        runs

        > successfully but produces an entirely black raster map as output

        > (Screen_02). All the pixels of this raster contain a value of 1.

        >

        > There is obviously something wrong here. Any tips? Should I use

        > “v.surf.idw” instead?

        >

        I found on an old maillist post that r.surf.idw works only
        with integer values. That's why you're getting everything
        "rounded" to 1. To work around this, you can try two options:

        * revert back to vectors, as you suggest. (v.surf.idw accepts
        decimal

        values)

             r.to.vect -z in=risk_raster out=risk_vector feature=point

             v.surf.idw in=risk_vector out=risk_idw layer=0

        OR

        * multiply your risk raster by some constant, say 1000, to get
        integer

        values. Then run r.surf.idw on that, and divide the result by
        the same

        1000 again to get back to your normalized values.

             r.mapcalc risk_integer=risk_raster*1000

             r.surf.idw in=risk_integer out=risk_idw_tmp

             r.mapcalc risk_idw=risk_idw_tmp/1000.0

             # Use the decimal value 1000.0 to insure that the result is

        floating point.

        Regards,

        Micha

        > Thank you so very much for your time Micha.

        >

        > Luciano

        >

        > *De:*Micha Silver [mailto:micha@arava.co.il
        <mailto:micha@arava.co.il>
        <mailto:micha@arava.co.il%20%3Cmailto:micha@arava.co.il%3E>]

        > *Enviado el:* Monday, November 05, 2012 4:01 AM

        > *Para:* Luciano La Sala

        > *Asunto:* [Bulk] Re: Two simple questions

        >

        > On 05/11/2012 04:42, Luciano La Sala wrote:

        >

        > Hello Micha,

        >

        > Where do you access the

        > modules///v.surf.idw/and///r.surf.idw/from?I just can’t
        seem to

        > find them!

        >

        > If you run GRASS on its own (not the QGIS plugin), these
        modules are

        > under the "Raster->Interpolate Surfaces"menu. In addition,
        all GRASS

        > commands can be run straight from the command line.

        >

        > Thanks!

        >

        > Luciano

        >

        >

        >

        > This mail was received via Mail-SeCure System.

        >

        >

        >

        >

        > --

        > Micha Silver

        > GIS Consulting

        > 052-3665918

        >http://www.surfaces.co.il

        >

        >

        > This mail was received via Mail-SeCure System.

        --

        Micha Silver

        GIS Consulting

        052-3665918

        http://www.surfaces.co.il

        This mail was received via Mail-SeCure System.

    This mail was received via Mail-SeCure System.

--
Micha Silver
GIS Consulting
052-3665918
http://www.surfaces.co.il

This mail was received via Mail-SeCure System.

--
Micha Silver
GIS Consultant, Arava Development Co.
http://www.surfaces.co.il