[raster] st_circle

[nclay]

A map algebra function that returns statistics on a circle defined by a radius.

[pracine]

What about ST_SummaryStats(ST_Clip(ST_Buffer()))?

[nclay]

This is a function to allow for statistics on different window shapes within map algebra instead of a square or rectangular window. See: ​http://resources.arcgis.com/en/help/main/10.1/index.html#//009z000000r7000000 My description was not all that clear.

Thanks,

Nathan

[Bborie Park]

I guessed this, st_wedge and st_annulus were related to doing focal stats. I wonder if this is better done as separate functions as currently ticketed or something more general. Essentially, these functions are masks. For my work I have a mapalgebra callback function that supersedes itself between mapalgebra and the actual processing callback function since I need to clip the neighborhood processed to a circle of some radius.

[nclay]

I've all ready written a function but have yet to test, that takes the map algebra window in and returns only the cells that fall within a r1,R2, thetaStart, thetaEnd:

 _st_focal( value double precision[][][],z integer, r1 integer, R2 integer, thetaStart double precision, thetaEnd double precision)

I ticked these functions as separate functions, we could use st_focal for all of them as they are really subsets of each other.

Yes, it could be implemented as a mask. Then the focal function would be something like:

 _st_focal( value double precision[][][],mask double precision[][][] )

It may be faster?

Thanks,

Nathan

[Bborie Park]

I personally prefer the flexibility of mask as that allows arbitrary areas of interest, e.g. ellipses or stars or diamonds.

[nclay]

So, we will need a function that returns a Mask of appropriate resolution for use in _st_focal to create convenience functions such as st_circle or st_annulus. Something like:

double precision [][][] ST_Mask(ST_CLIP(RASTER,GEOMETRY,BITMAP or WEIGHT))

That way we could also implement weighted masks, along with bit map masks.

Thanks,

Nathan

[Bborie Park]

Yep. Something like that. See if it works…

[nclay]

I have already coded up a mask function that takes in a geometry and a x, y resolution and produces a mask. That is all well and good if you know what resolution that you want you mask to be. It would be nice to be able to pass in a geometry and a SRID and have it pass out an appropriately sized mask. This would enable, one to pass in a Geometry describing distances eg a circle 16 meters wide and have it project on to an appropriate resolution mask for that SRID.

Can you think of anyway of doing this currently?

[Bborie Park]

If the geometry is passed along with a raster of the same SRID, the spatial grid of the raster can be used to convert the geometry to a mask. Without that, there is no way to know how to create a correct mask.

So…

WITH foo AS (
	SELECT ST_Mask(ST_AsRaster(geomtomask, refraster)) AS mask
	FROM myrasttable
	LIMIT 1
)
SELECT
	...
FROM foo

You could call ST_AsRaster() on the geometry without any reference raster but you'd have to align the output raster of ST_AsRaster() with ST_Resample() before you could actually use it. But assuming that a mask is being generated to use against an existing set of rasters, this shouldn't be an issue.

[nclay]

Mangled VARIADIC? I am calling map algebra as follows:

RETURN ST_MapAlgebra(
			ARRAY[ROW(_rast, _nband)]::rastbandarg[],
			'_st_focal4ma(double precision[][][], integer[][], text[])'::regprocedure,
			_pixtype,
			_extenttype, _customextent,
			r1, r1,VARIADIC ARRAY[ROW(_mask,op)]::text[]);

Yet inside st_focal4ma I am getting the following for userargs:

userarg  {"(\"{{{NULL,NULL,NULL,1,1,1,1,NULL,NULL,NULL},{NULL,1,1,1,1,1,1,1,1,NULL},{NULL,1,1,1,1,1,1,1,1,NULL},{1,1,1,1,1,1,1,1,1,1},{1,1,1,1,1,1,1,1,1,1},{1,1,1,1,1,1,1,1,1,1},{1,1,1,1,1,1,1,1,1,1},{NULL,1,1,1,1,1,1,1,1,NULL},{NULL,1,1,1,1,1,1,1,1,NULL},{NULL,NULL,NULL,1,1,1,1,NULL,NULL,NULL}}}\",var)"}

This is not a valid array, is this of my doing, a bug or expected?

Thanks,

Nathan

[Bborie Park]

userarg is a valid array with 2 text elements.

Your VARIADIC ARRAY[…]::text[] is no different than…

RETURN ST_MapAlgebra(
			ARRAY[ROW(_rast, _nband)]::rastbandarg[],
			'_st_focal4ma(double precision[][][], integer[][], text[])'::regprocedure,
			_pixtype,
			_extenttype, _customextent,
			r1, r1, _mask::text, op::text);

[nclay]

I've implemented st_circle in plpgsql and I am not impressed with its speed.

So I am considering modifying nMapAlgebra to accept a mask argument. A mask would have the following properties: height,width, and weighted. If weighted is marked the mask is used to weight the pixel values as: pixval * mask. If it is not marked then the pixel values are returned where the mask is not null.

To this end, I think it would be best to modify nMapAlgebra by replacing rt_band_get_pixel with rt_band_mask_get_pixel this function would consider the mask when returning pixel values. Do you (dustymugs) have any suggestions or objections? Should I open another ticket for this?

struct rt_mask_t {
  uint16_t width;
  uint16_t height;
  bool weighted;
  void** data;
};

/**
 * Get Pixel from band with consideration of a masked extent returning nodata value
 * if mask position is null.
 *@param band : the band to get the pixel value from
 *@param mask : the mask to mask the band with
 *@param extent : array of values containing the bounds of the mask.
 *@param x : the x coordinate of the band to be sampled
 *@param y : the y coordinate of the band to be sampled
 *@param pixval : the value of the band to be returned
 *@param nodata : flag to indicate that the value is nodata or masked
 *
 *@return ES_ERROR on error
 */
int rt_band_mask_get_pixel(rt_band band, 
			   rt_mask mask,
			   int** extent,
			   int x,
			   int y,
			   double *pixval,
			   int* nodata
) 

Thanks,

Nathan

[Bborie Park]

I'm not surprised you went this way for performance. The rt_band_mask_get_pixel should really get all the necessary pixel values at once and be a wrapper around rt_band_get_pixel and/or rt_band_get_pixel_line.

rt_errorstate rt_band_get_pixels(
  rt_band band
  rt_mask mask,
  int x,
  int y,
  double **pixel,
  int **nodata
)

I don't get the purpose of extent since rt_mask should have the requisite information to describe itself.

Also, how do you envision the function call at the SQL level? I ask because you'll want to create a ticket for that function call to which we can move this discussion.

[nclay]

dependant on st_wedge private function _st_wedge

[robe]

Milestone renamed