I am working on a problem where I have polygons for features and I want look at daily data from MODIS on, and I am looking for suggestions on how to determine which pixels it is that I need to look at. I have tried gdallocationinfo however this only gives me a single pixel, and I am looking to get all the pixels that are at least 80% inside of the polygon. I am looking for a coding based solution as I have about 8,000 polygons I'm looking to process.

One other thing that makes the process a bit more complicated, the MODIS bands I'm using at at 500m resolution, while the cloud-mask for it is at 1km (int the state_1km band) so I will need to keep the two data-sets in alignment throughout the process.


The solution I came up with was to put one of the rasters into postGIS with raster2pgsql, and then use this script, which took about 9 hours to run.

INSERT INTO modis.pixels (SELECT t1.v, t2.v, ST_PixelASPolygon(raster.rast, t1.v, t2.v)
FROM modis.raster
CROSS JOIN generate_series(1, 2400) as t1(v)
CROSS JOIN generate_series(1, 2400) as t2(v));

Then I built an index which took about 2 and a half minutes, and finally used this script which was taking about 30ms/polygon.

INSERT INTO modis.polygon_pixels(
SELECT polygon_id::integer, 'h10v03'::text AS tile, x::smallint, y::smallint,(ST_Area(ST_TRANSFORM(ST_Intersection(st_transform(boundary,96842), geom), 3857))/ST_AREA(ST_TRANSFORM(geom, 3857)) * 100)::smallint AS coverage
FROM polygon.boundary
JOIN modis.pixels
    ON ST_Intersects(st_transform(boundary,96842), geom));

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.