Extracting exact are of a pixel being overlapped by a JimVec object

[kasra-keshavarz]

I am wondering if there is a way to extract the exact area of each pixel (of a GeoTIFF file loaded as a Jim object) being overlapped by a JimVec object (i.e., a polygon ESRI shapefile).

[kempenep]

Hi Kasra.

I don't think it can easily be done with a Jim raster object directly. Perhaps there is a different end-goal behind your question that can be handled differently. But then I would have to understand the context better. The problem is that the vector is first rasterized with the extract functions in pyjeo. One way might be to up-sample the raster, but this is not very efficient and still an approximation.

Perhaps a better way might be to vectorize the pixels and calculate the intersection with the vector. This can be done as suggested here

from osgeo import ogr

wkt1 = "POLYGON ((1208064.271243039 624154.6783778917, 1208064.271243039 601260.9785661874, 1231345.9998651114 601260.9785661874, 1231345.9998651114 624154.6783778917, 1208064.271243039 624154.6783778917))"
wkt2 = "POLYGON ((1199915.6662253144 633079.3410163528, 1199915.6662253144 614453.958118695, 1219317.1067437078 614453.958118695, 1219317.1067437078 633079.3410163528, 1199915.6662253144 633079.3410163528)))"

poly1 = ogr.CreateGeometryFromWkt(wkt1)
poly2 = ogr.CreateGeometryFromWkt(wkt2)

intersection = poly1.Intersection(poly2)

print intersection.ExportToWkt()

But then it also depend on how many pixels you would like to extract.

[kasra-keshavarz]

Thanks for your input. The idea is accurate zonal statistics for the purpose of hydrological modelling on a large scale if you are interested in the end goal. In this case, the number of pixels can grow exponentially with an increase in scale, so I was looking for an efficient solution to do so.

I know exactextract is offering a similar, efficient solution, but the jeolib is accessible from Python and would've been a much flexible tool if it offers such a feature.