Open kasra-keshavarz opened 2 years ago
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.
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.
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).