[GRASS-user] Still in trouble
Micha Silver
micha at arava.co.il
Sat Nov 10 07:38:42 PST 2012
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 at 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 at 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 at 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 at arava.co.il
> <mailto:micha at arava.co.il>
> <mailto:micha at arava.co.il%20%3Cmailto:micha at 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
More information about the grass-user
mailing list