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 <emiliomayo...@gmail.com> wrote: > n Wed, Aug 26, 2009 at 9:37 AM, Jose Gomez-Dans<jgomezd...@gmail.com> wrote: >> 2009/8/26 Dylan Beaudette <debeaude...@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 <debeaude...@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@lists.osgeo.org > http://lists.osgeo.org/mailman/listinfo/gdal-dev > _______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
