| Re: How to extract pixel values from a GeoTIFF using an Esri Shapefile [message #91745 is a reply to message #74179] |
Wed, 19 August 2015 14:57   |
Adam Erickson
Messages: 8
Registered: July 2015
|
Junior Member |
|
|
On Monday, January 3, 2011 at 4:25:55 PM UTC-8, David Fanning wrote:
> Guillermo writes:
>
>> ; populate the mask (assumes is in the same projection
>> ;and covers the same extent as the geotiff)
>> myshape->IDLffShape::GetProperty, N_ENTITIES=n
>> FOR i=0L, n-1 DO BEGIN
>> feati= myshape->IDLffShape::GetEntity(i)
>> featix= Round((Reform((*feati.vertices)[0,*])-x0)/psz)
>> featiy= Round((y0-Reform((*feati.vertices)[1,*]))/psz)
>> featis= POLYFILLV(featix, featiy, ns, nl)
>> IF featis[0] NE -1 THEN mask[featis]= feati.ishape +1
>> ENDFOR
>
> This is nearly identical to the solution I sent Paul
> earlier today. The problem with IDLanROI is that it
> can't be used in the native projected meter coordinates
> that the image and shape files are in. The projected
> meter coordinates have to be "converted" to "pixel"
> coordinates by subtracting the offset and dividing
> by the image range, before they can be loaded in the
> object.
>
> I didn't round my values, and I don't think you need
> to do so here, even with PolyFillV. In fact, I think
> you might get slightly more accurate values by not
> rounding, although this is a quibble with Paul's
> image.
>
> I've made myself a note to write an article when
> I get some time. :-)
>
> Cheers,
>
> David
>
>
> --
> David Fanning, Ph.D.
> Fanning Software Consulting, Inc.
> Coyote's Guide to IDL Programming: http://www.dfanning.com/
> Sepore ma de ni thui. ("Perhaps thou speakest truth.")
Honestly, use QGIS for this. It's 64-bit, unlike ArcGIS, and its GDAL functions outperform everything else I've tried, aside from my own GDAL python code. Simply add the Zonal Statistics tool in Plugins -> Manage and Install Plugins... then select Zonal statistics under the Raster menu. The operation computes in one minute on very large mosaics with many polygons and writes the values to new fields in the shapefile attribute table - something I could not get ENVI to do.
ArcGIS ran for three days before crashing (one bad polygon is all it took), with the raster placed in a File Geodatabase to mitigate its 32-bit limitations. ENVI took a day to compute this, plus one hour simply to draw the polygons to screen, while outputting the data in what appears to be a useless format. Originally, I was doing this operation with IDL similar to the above methods, but my mosaic is much larger than my 96 GB of RAM. My web browser uses more RAM than QGIS, which draws to screen very quickly as well. I cannot recommend it enough for zonal statistics.
Cheers,
Adam Erickson
Integrated Remote Sensing Studio
UBC
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
