Wednesday, September 13, 2017

Compute zonal statistics on area of interest polygon on fly

Draft to compute zonal statistics on area of interest polygon or shape draw on leaflet using PostgreSql and Java -

Geometry extracted from leaflet as GeoJson.

Polygon-
{"type":"Feature","properties":{"desc":null,"image":null},"geometry":{"type":"Polygon","coordinates":[[[-117.103271484375,34.264026473152875],[-117.1142578125,34.14818102254435],[-117.03186035156251,34.10498222546687],[-116.91925048828124,34.14363482031264],[-116.94946289062499,34.25494631082515],[-117.0867919921875,34.252676117101515],[-117.103271484375,34.264026473152875]]]}}


{"type":"Feature","properties":{"desc":null,"image":null},"geometry":{"type":"Polygon","coordinates":[[[-116.86981201171875,34.11407854333859],[-116.86981201171875,34.24132422972854],[-116.71874999999999,34.24132422972854],[-116.71874999999999,34.11407854333859],[-116.86981201171875,34.11407854333859]]]}}

Point-
{"type":"Feature","properties":{"desc":null,"image":null},"geometry":{"type":"Point","coordinates":[-116.9879150390625,34.048108084909835]}}


Area of GeoJson
Select (ST_Area(ST_GeomFromText
('POLYGON ((-117.16918945312501 34.27083595165,-117.1307373046875 34.166363384737892,-117.00988769531251 34.161818161230386,-116.93298339843751 34.239053668516412,-116.98516845703126 34.27083595165,-117.16918945312501 34.27083595165))',4326)))



Zonal statistics for each raster tile -

SELECT rid, (ST_SummaryStats (ST_Clip(rast,ST_GeomFromText('POLYGON ((-117.16918945312501 34.27083595165,-117.1307373046875 34.166363384737892,-117.00988769531251 34.161818161230386,-116.93298339843751 34.239053668516412,-116.98516845703126 34.27083595165,-117.16918945312501 34.27083595165))',4326),true))).*
FROM public.biodiv_ssolnw_wgs84
WHERE ST_Intersects
(rast,ST_GeomFromText
('POLYGON ((-117.16918945312501 34.27083595165,-117.1307373046875 34.166363384737892,-117.00988769531251 34.161818161230386,-116.93298339843751 34.239053668516412,-116.98516845703126 34.27083595165,-117.16918945312501 34.27083595165))',4326))
Group By rid



Combined zonal statistics for raster tiles as a single raster -

SELECT  (ST_SummaryStats (ST_Clip(ST_Union(rast), ST_GeomFromText('POLYGON ((-117.16918945312501 34.27083595165,-117.1307373046875 34.166363384737892,-117.00988769531251 34.161818161230386,-116.93298339843751 34.239053668516412,-116.98516845703126 34.27083595165,-117.16918945312501 34.27083595165))',4326),true))).*
FROM public.biodiv_ssolnw_wgs84
WHERE ST_Intersects
(rast,ST_GeomFromText
('POLYGON ((-117.16918945312501 34.27083595165,-117.1307373046875 34.166363384737892,-117.00988769531251 34.161818161230386,-116.93298339843751 34.239053668516412,-116.98516845703126 34.27083595165,-117.16918945312501 34.27083595165))',4326))


Function to Compute Zonal Statistics -
private String GetSqlQuery(String tableName,String polygonGeom){
        String pgSqlQuery =  "SELECT row_to_json(t) from ( SELECT (ST_SummaryStats (ST_Clip(ST_Union(rast), ST_GeomFromText('"+polygonGeom+"',4326),true))).*,"
                +" ST_Area(ST_GeomFromText('"+polygonGeom+"',4326))"
                +" FROM public."+tableName
                +" WHERE ST_Intersects(rast,ST_GeomFromText ('"+polygonGeom+"',4326)))t";

        return pgSqlQuery;
    }

I will make full post soon.

/* polygon 1 */
SELECT  (ST_SummaryStats (ST_Clip(ST_Union(rast), ST_GeomFromText('MULTIPOLYGON (((-117.603149 34.034453,-117.603149 34.079962,-117.526245 34.079962,-117.526245 34.034453,-117.603149 34.034453)))',4326),true))).*
FROM public.biodiv_ssolnw
WHERE ST_Intersects(rast,ST_GeomFromText('MULTIPOLYGON (((-117.603149 34.034453,-117.603149 34.079962,-117.526245 34.079962,-117.526245 34.034453,-117.603149 34.034453)))'  ,4326));

/* polygon 2 */
SELECT  (ST_SummaryStats (ST_Clip(ST_Union(rast), ST_GeomFromText('MULTIPOLYGON (((-119.274445 34.644507,-119.274445 34.663711,-119.251099 34.663711,-119.251099 34.644507,-119.274445 34.644507)))',4326),true))).*
FROM public.biodiv_ssolnw
WHERE ST_Intersects(rast,ST_GeomFromText('MULTIPOLYGON (((-119.274445 34.644507,-119.274445 34.663711,-119.251099 34.663711,-119.251099 34.644507,-119.274445 34.644507)))'  ,4326));

/* both */
SELECT  (ST_SummaryStats (ST_Clip(ST_Union(rast), ST_GeomFromText('MULTIPOLYGON (((-119.274445 34.644507,-119.274445 34.663711,-119.251099 34.663711,-119.251099 34.644507,-119.274445 34.644507)),((-117.603149 34.034453,-117.603149 34.079962,-117.526245 34.079962,-117.526245 34.034453,-117.603149 34.034453)))',4326),true))).*
FROM public.biodiv_ssolnw

WHERE ST_Intersects(rast,ST_GeomFromText('MULTIPOLYGON (((-119.274445 34.644507,-119.274445 34.663711,-119.251099 34.663711,-119.251099 34.644507,-119.274445 34.644507)),((-117.603149 34.034453,-117.603149 34.079962,-117.526245 34.079962,-117.526245 34.034453,-117.603149 34.034453)))'  ,4326));

Leaflet , PostGIS , Postgres

0 comments :

Post a Comment

 

© 2011 GIS and Remote Sensing Tools, Tips and more .. ToS | Privacy Policy | Sitemap

About Me