[postgis-users] Raster pixel count too high

Sun Jun 23 12:52:59 PDT 2013

Hello,

I have a lot of work and I won't be able to test your data now, but I'll try ASAP.

On my side, I have to deal with rasters dem tiles which overlap each other. I need to have unique tiles and I made the two attached functions to achieve that :
1. The first one, makegrid, simply builds a grid given a bounding box geometry and height / width
2. The second, unique_tiles, calls the first one to make a whole raster made of regular tiles of the wanted size in number of pixels from a raster table of aligned tiles (it finds its bounding box and use it to makegrid from left bottom corner using original pixel height and width and then use this grid to clip raster)

Maybe they can help you in your clipping process. You may improve performance if you remove the st_union in unique_tiles (needed if your raster is made of multiple tiles).

Hugues.

De : postgis-users-bounces at lists.osgeo.org [mailto:postgis-users-bounces at lists.osgeo.org] De la part de Andreas Forø Tollefsen
Envoyé : dimanche 23 juin 2013 10:47
À : PostGIS Users Discussion
Objet : Re: [postgis-users] Raster pixel count too high

Thanks for your answers.

First to Hugues. I do not think they are perfectly aligned. The raster we have imported seems to start just west of -180. Hence, it is not within the limits of SRID 4326. We have to modify it a little.

What I do find strange is that ArcGIS counts 3600 which is the expected pixel count. So why these two functions count so differently on the exact same data is weird.

Try yourself with one of the Nightlights data: http://www.ngdc.noaa.gov/eog/dmsp/downloadV4composites.html <http://www.ngdc.noaa.gov/eog/dmsp/downloadV4composites.html>  and our vector grid shapefile available here: http://www.prio.no/Data/PRIO-GRID/ <http://www.prio.no/Data/PRIO-GRID/> 

Kim: We do not directly use ST_Intersection in our script, but ST_Clip to clip the raster according to our polygons, and only where they ST_Intersects(), not ST_Intersection(); 

I still do not understand your suggestion.

Best,

Andreas

2013/6/18 Kim Bisgaard <kib at dmi.dk <mailto:kib at dmi.dk> >

Hi,

Because ST_Intersection() returns neighbour pixels sharing a same value as only one
 polygon. You thus get a bigger area and thus more points.

Bit me once :-/

Regards,
Kim

On 2013-06-18 14:59, Andreas Forø Tollefsen wrote:

	Hi Kim, 

	Thanks for your answer. However, we want raster as an output, since we want to be able to use the summarystats function.

	Please elaborate how you think ST_PixelAsPolygons should solve out issue?

	Thanks.

	Andreas

	2013/6/18 Kim Bisgaard <kib at dmi.dk <mailto:kib at dmi.dk> >

	Hi,

	Try to use 'ST_PixelAsPolygons(ST_Clip(n.rast, p.cell))'
	instead of 'ST_Intersects(n.rast, p.cell)'

	Regards,
	Kim 

	On 2013-06-18 11:03, Andreas Forø Tollefsen wrote:

		Hi, 

		We are working on a raster summarystats script to calculate various statistics for the pixels within fishnet polygons.

		Our raster cell size is 0.0083333333333... x 0.0083333333333... degrees while our quadrat polygons are 0.5 x 0.5 decimal degrees.

		This should give us 60x60 raster pixels within each of our polygons. ArcGIS zonal statistics returns a pixel count of 3600 in addition to other statistics.

		However, PostGIS returns 3721 pixel count.

		We do not really understand why, but it seems that our query includes some pixels that are outside of the polygon, but still touches the vertices of the polygon and are therefore included in the calculation.

		Are there any way of modifying our script to return the same result as ArcGIS?

		Thanks!

		Andreas

		script:

		/* This query makes one raster for each PRIO-GRID cell. Clip and union is the procedure. */

		INSERT INTO nightlightsprio (gid, "year", rast)

		(SELECT gid, "year", ST_Union(raster) as rast

		FROM 

		(SELECT p.gid, n."year", ST_Clip(n.rast, p.cell) as raster 

		FROM nightlights n, priogridyear p

		WHERE ST_Intersects(n.rast, p.cell)

		AND n."year" = p."year"

		) 

		as priorast

		GROUP BY gid, "year");

		/* Default BandNoDataValue is 0. Raster value 0 means no light, not no data. Setting to NULL. This produces correct results. */

		UPDATE nightlightsprio2 SET rast = ST_SetBandNoDataValue(rast, 1, NULL);

		ALTER TABLE nightlightsprio2 ADD COLUMN nightlights_sum double precision,

		ADD COLUMN nightlights_mean double precision,

		ADD COLUMN nightlights_sd double precision,

		ADD COLUMN nightlights_min double precision, 

		ADD COLUMN nightlights_max double precision, 

		ADD COLUMN nightlights_count integer;

		UPDATE nightlightsprio2 SET nightlights_sum = (ST_SummaryStats(rast)).sum;

		UPDATE nightlightsprio2 SET nightlights_mean = (ST_SummaryStats(rast)).mean;

		UPDATE nightlightsprio2 SET nightlights_sd = (ST_SummaryStats(rast)).stddev;

		UPDATE nightlightsprio2 SET nightlights_min = (ST_SummaryStats(rast)).min;

		UPDATE nightlightsprio2 SET nightlights_max = (ST_SummaryStats(rast)).max;

		UPDATE nightlightsprio2 SET nightlights_count = (ST_SummaryStats(rast)).count;

		_______________________________________________
		postgis-users mailing list
		postgis-users at lists.osgeo.org <mailto:postgis-users at lists.osgeo.org> 
		http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users <http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users> 

	-- 
	Kim Bisgaard

	Application Development Division     Phone: +45 3915 7562 <tel:%2B45%203915%207562>  (direct)
	Danish Meteorological Institute      Fax: +45 3915 7460 <tel:%2B45%203915%207460>  (division)

	_______________________________________________
	postgis-users mailing list
	postgis-users at lists.osgeo.org <mailto:postgis-users at lists.osgeo.org> 
	http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users <http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users> 

	_______________________________________________
	postgis-users mailing list
	postgis-users at lists.osgeo.org <mailto:postgis-users at lists.osgeo.org> 
	http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users <http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users> 

-- 
Kim Bisgaard

Application Development Division     Phone: +45 3915 7562 <tel:%2B45%203915%207562>  (direct)
Danish Meteorological Institute      Fax: +45 3915 7460 <tel:%2B45%203915%207460>  (division)

_______________________________________________
postgis-users mailing list
postgis-users at lists.osgeo.org <mailto:postgis-users at lists.osgeo.org> 
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users <http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users> 

-------------- next part --------------
A non-text attachment was scrubbed...
Name: function_makegrid.sql
Type: application/octet-stream
Size: 1282 bytes
Desc: function_makegrid.sql
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20130623/d7bb420a/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: function_unique_tiles.sql
Type: application/octet-stream
Size: 3875 bytes
Desc: function_unique_tiles.sql
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20130623/d7bb420a/attachment-0001.obj>