[gdal-dev] calculate raster values within a vector

[Roger André]

Just to add my 2 cents on this.  I've done this sort of thing in the
past by using gdal_rasterize to set all pixels which are outside the
polygon(s) of interest to nodata. Then I use a very simplistic Python
script to dump all of the pixel values to STDOUT and redirect to a csv
file.  It would be pretty simple to trap the values instead of dumping
them, and do any sort of stat query you want in Python.  I was shown
an excellent resource for learning to use gdal to do this type of work
- http://www.gis.usu.edu/~chrisg/python/2009/.

I believe that if I had to do it now, I would probably just load them
all up into GRASS, run v.to.rast on the shapefile to create a
rasterized version of the polygons, and then use masks, or r.mapcalc,
in conjunction with r.univar to generate the statistics.

I totally hear what you mean though, Emilio. It's nice when your
coworkers can access and use your work too.  I suspect that factor
will determine ultimately which approach is best suited for the task.
I do highly recommend Mitasova's book on GRASS though.  It will have
you up and running in no time.

Roger
--



On Wed, Aug 26, 2009 at 12:02 PM, Emilio Mayorga
<emiliomayorga at gmail.com> wrote:
> n Wed, Aug 26, 2009 at 9:37 AM, Jose Gomez-Dans<jgomezdans at gmail.com> wrote:
>> 2009/8/26 Dylan Beaudette <debeaudette at ucdavis.edu>
>>>
>>> ndeed. The current (C++) version of Starspan is more or less unsupported.
>>> I tend to use an older, stabler, version-- but have become frustrated with
>>> it in recent studies. I think that the current maintainer Jon Greenberg is
>>> working on an R implementation-- however I am not certain that this will
>>> scale well to very large rasters. A Python incantation would be more
>>> flexible
>>> than the current C++ version, but perhaps at a speed cost. The current
>>> version is blazing fast, but there aren't any C++ programmers working on
>>> it now. If there is sufficient interest, I would like to get it into OSGeo so
>>> that a more skilled programmer (than myself) can have a look at it.
>>
>> Python is extremely fast for some of these things. The version I tend to use
>> < http://sites.google.com/site/spatialpython/zonal-statistics> is bearable
>> even working with whole landsat scenes. The example there (with one MODIS
>> 1km tile) takes in my laptop:
>> %timeit S = ZonalStats ( labels_file, data_file, stats='mean') ; std_s =
>> S.ResampleToGrid("std")
>> 10 loops, best of 3: 1.44 s per loop
>>
>> %timeit S = ZonalStats ( labels_file, data_file, stats='mean') ; std_s =
>> S.ZonalStatistics("mean")
>> 10 loops, best of 3: 359 ms per loop
>>
>>
>> I plan to add a number of refinements to the code, such as adding a vector
>> input that gets automaticall rasterised, something that the GDAL python
>> wrappers have allowed for quite some time now, and adding these values to
>> the output as extra fields, but so far, I have had no need for these, so I
>> haven't really bothered.
>
> Wow, Jose, thanks for sharing that! It looks terrific. I may be able
> to contribute to it, but not until October or so.
>
>
>> 2009/8/26 Dylan Beaudette <debeaudette at ucdavis.edu>
>> If you are doing raster-on-raster zonal stats, then GRASS is fairly good at
>> that too.
>
> Dylan, thanks for your comments on Starspan. Regarding GRASS, I wish I
> had time to start using it, but so far I haven't been able to. One of
> these days... Still, accessing this functionality from Python is a big
> plus for me. For example, I'm sharing some of my zonalstats-like code
> with a colleague. I've helped him set up Python with all the right
> modules, and he can run my code. Adding GRASS to that bag of new
> things would probably make this task too complicated. Anyways, that's
> just my own preferred workflow. I'll get to GRASS one of these days...
>
> Cheers,
>
> -Emilio
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>