[GRASS-user] Re: count pixels by attribute in sampling units
Patrick Giraudoux
patrick.giraudoux at univ-fcomte.fr
Sat Dec 29 04:00:37 EST 2007
Hamish wrote:
> [r.le.setup]
>
> Hamish wrote:
>
>> bug found and fixed in 6.3 svn. The logic around the >100 check was
>> to blame. Thanks for pin-pointing it.
>>
>>
> http://trac.osgeo.org/grass/log/grass/trunk/raster/r.le/r.le.setup/sample.c
>
> Patrick:
>
>> OK I will do it asap. However I don't know how to insert the new code
>> practically from the script given
>>
> ....
>
>> However, I don't know how to insert the C code in a GRASS working
>> version (I can manage with a patch and a cvs version since recent
>> only...
>>
>
> You can get a patch from the above link. Click on the 'View Changes'
> button, then the 'Unified diff' link. Or from the changes page just see
> how to edit sample.c by hand.
>
>
>> (I am using R and grass but I am not a C developper)...
>>
>
> You don't have to know any C to apply the patch, nor much C to
> understand the change. Simply knowing if{;} else{;} syntax and how to
> /* comment out lines */ is enough to start tinkering with the code.
>
> In this case the "else" part was only run if the radius was under 100,
> instead of being run after radius of any size was accepted.
>
>
>> Actually I am a landscape ecologist and R programmer.
>>
>
> Well FWIW, I'm actually an engineer and oceanographer in a place
> dominated by smart marine biologists and ecologists. I think only a
> couple of grass developers are real programmers, and their ongoing
> patience with the rest continues to amaze me.
>
>
>> Has it been inserted in the last cvs version (guess not from your
>> mail) at http://grass.itc.it/grass63/source/snapshot/, or any other
>> way ?
>>
>
> The fix has been applied in 6.3-svn, not CVS. In the last weeks GRASS
> has migrated from CVS to a new SVN source code repository at OSgeo.
> Sorry the source code snapshot links on the website are currently out
> of date, check a recent post by Markus on the grass-dev mailing list
> for a SVN-snapshot link. I am done trying to do source code commits
> from the road, so someone else will have to update that or wait for me
> to get back to the office.
>
> Here is a link to the unidiff if you don't mind patches:
> http://trac.osgeo.org/grass/changeset?format=diff&new=29503&old=29458&new_path=grass%2Ftrunk%2Fraster%2Fr.le%2Fr.le.setup%2Fsample.c&old_path=grass%2Ftrunk%2Fraster%2Fr.le%2Fr.le.setup%2Fsample.c
>
> Patrick:
>
>> Sorry for your virus mishappenings...
>>
> ....
>
>> but took 15 days to solve the trouble.
>>
>
> I guess my lesson is that there is no such thing as a "safe" MS Windows
> machine, at least with any confidence. It wasn't a major infection or
> anything, just some internet-explorer auto-installed spyware. It was
> easily found and removed, but enough to spook me.
>
>
>
>> 2/ Regarding r.le*, it would be nice to get simple and quick reports
>> on sampling unit composition (e.g. number of pixels of each attribute
>>
>
>
>> category, mean, min, max, etc. e.g. the same as for r.neighbors)
>> without going straight to time-consuming shape analysis. Indeed,
>> classical neighborhood analysis is limited to 25 pixels in GRASS and
>> there are many cases where statistics on larger sampling units are
>> needed. It may be that such simple statistics on sampling units can
>> be obtained with other GRASS functions, but I cannot find one
>> throughout books and threads on the web... and no feed back yet on
>> this on the GRASS list.
>>
>
> it is an easy edit of the r.neighbors source code if you want to make
> that number more than 25.
>
> but yes, there are modules to do this.
> r.buffer + r.univar,
> v.buffer + v.rast.stats
> v.what.rast
> .... more?
>
> It might be nice for that to be built into r.le.* automatically, and
> probably not that hard to do, but generally I would prefer if new
> development went into r.li, and r.le only got bug fixes. I don't mind
> simple niceties being added but it seems like wasted effort when r.li
> is supposed to replace r.le. Of course it is open source so if someone
> prefers to keep adding to r.le I'm not going to stop them or refuse
> patches.
>
>
>> So, no real trouble with the doc, which looks fair enough.
>>
>
> In general for bugs I would rather fix them than document them.
> (but document the "it's a feature not bug" design choices)
>
>
>
> regards,
> Hamish
>
>
Dear Hamish,
I have modified the code by hand (no patch used, indeed this was quite
easy...) and recompiled a cvs version. Everything works fine now in
r.le.setup with sample unit radius > 100. This gave me the opportunity
to see how the file system is organised in GRASS and may encourage me to
go further... Thanks for the incent !
One question remains pending about r.li. It does not accept any circle
radius at the moment for some reasons (see the threat 'r.li.setup
trouble with circle').
Regarding crude totals in sampling unit circles, I have written a
UNIX/GRASS 'pilot' programme which does what I need. I though it may
help other listers and thus you'll find the script below.
Cheers,
Patrick
#!/bin/sh
##########################################################################
# Program to run within GRASS (this file must be copied into a folder #
# of the PATH). A circular buffer is created from each point location #
# of a vector file and it creates a text file of two columns: #
# col 1 = pixel value, col 2 = pixel count corresponding to this value. #
# Buffer results are separated by '$ $'. This file can be read from #
# e.g. R and formatted for user purposes. #
# #
# The programme shall be used on isometric projected coordinates #
# (UTM, Lambert, etc...) - the size of the buffer does not adapt to #
# geographical coordinates in degrees. #
# #
# Before running the programme, do not forget: #
# - to name the vector file (variable fichier=filename1, see below) #
# - to name the raster file (variable rast=filename2, see below) #
# - to give the buffer radius (in kilometres) (variable radius=) #
# #
# By defaut results are found in a file whose name is the name #
# of the point file + the suffix .txt #
# #
# This program is just a prototype and few controls are really done. For #
# instance, users must make sure that point vector file and raster file #
# do overlap; null values are however counted. #
# #
# P. Giraudoux - version 28.12.2007 21:11:40 #
##########################################################################
if test "$GISBASE" = ""; then
echo "You must be in GRASS to run this programme."
exit
fi
fichier=chinaSP4 # name of the point file (vect)
rast=dem3 # name of the raster file (rast)
radius=100 # buffer radius (kilometers)
resultat=$fichier.txt # name of the result file
cote=`expr $radius \* 1000` # half of the length of the square edge of
the GRASS
# region (in meters) where the buffers will
take place
rm $resultat # remove an older file of the same name if any
points=`v.info -t $fichier | grep 'points=' | cut -d= -f2` # number of
point locations
count=0
while [ $count -lt $points ]
do
count=`expr $count + 1`
echo;echo
echo $count/$points
echo;echo
g.remove MASK --q # remove mask if any (some troubles happened with
the new mask when not done)
v.extract -t in=$fichier out=essais list=$count type=point --o --q #
point extraction
# get point coordinates and setup the GRASS region at buffer size
lat=`v.info -g essais | grep 'north=' | cut -d= -f2`
long=`v.info -g essais | grep 'east=' | cut -d= -f2`
nn=`echo "scale=2; $lat + $cote" | bc` # use bc for floating point
computing
ns=`echo "scale=2; $lat - $cote" | bc`
ne=`echo "scale=2; $long + $cote" | bc`
nw=`echo "scale=2; $long - $cote" | bc`
g.region n=$nn s=$ns e=$ne w=$nw # setup GRASS region
#-----------------------------------------
v.to.rast in=essais out=essais use=val type=point layer=1 value=1
--overwrite --q # convert the single point file into a raster
r.buffer in=essais out=essais.buf distances=$radius units=kilometers
--overwrite --q # computes the buffer
r.mask -o in=essais.buf maskcats=* --q # mask excluding the area not
in the buffer
#d.erase
#d.rast dem3
r.stats -c $rast >> $resultat # adds the result at the end of the
result file
echo "$ $" >> $resultat # buffer separator
done
g.remove MASK --q
g.remove rast=essais --q
g.remove rast=essais.buf --q
g.remove vect=essais --q
g.region rast=$rast
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20071229/2d74a7b2/attachment-0001.html
More information about the grass-user
mailing list