Opened 7 years ago

# [raster] st_circle

Reported by: Owned by: nclay Bborie Park medium PostGIS Fund Me raster master

### Description

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

### comment:2 Changed 7 years ago by 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

### comment:3 Changed 7 years ago by 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.

### comment:4 Changed 7 years ago by 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

### comment:5 Changed 7 years ago by Bborie Park

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

### comment:6 Changed 7 years ago by 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

### comment:7 Changed 7 years ago by Bborie Park

Yep. Something like that. See if it works...

### comment:8 Changed 7 years ago by 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?

### comment:9 Changed 7 years ago by 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 (
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.

### comment:10 Changed 7 years ago by 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,
```

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

### comment:11 Changed 7 years ago by 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);
```

### comment:12 Changed 7 years ago by 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** extent,
int x,
int y,
double *pixval,
int* nodata
)
```

Thanks,

Nathan

### comment:13 Changed 7 years ago by 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
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.