PostGIS  2.2.8dev-r@@SVN_REVISION@@

◆ RASTER_nMapAlgebra()

Datum RASTER_nMapAlgebra ( PG_FUNCTION_ARGS  )

Definition at line 508 of file rtpg_mapalgebra.c.

References ovdump::band, rtpg_nmapalgebra_arg_t::callback, rtpg_nmapalgebra_arg_t::cextent, rt_mask_t::dimx, rt_mask_t::dimy, rtpg_nmapalgebra_arg_t::distance, ES_NONE, ET_CUSTOM, ET_LAST, ET_SECOND, rtpg_nmapalgebra_arg_t::extenttype, FALSE, rtpg_nmapalgebra_arg_t::hasband, rtpg_nmapalgebra_arg_t::hasnodata, rtpg_nmapalgebra_arg_t::mask, rtpg_nmapalgebra_arg_t::nband, rt_iterator_t::nband, rt_iterator_t::nbnodata, rt_mask_t::nodata, rtpg_nmapalgebra_arg_t::nodataval, rtpg_nmapalgebra_arg_t::numraster, rtpg_nmapalgebra_arg_t::pgcextent, rtpg_nmapalgebra_arg_t::pixtype, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, PT_END, rtpg_nmapalgebra_arg_t::raster, rtrowdump::raster, rt_iterator_t::raster, rt_band_get_hasnodata_flag(), rt_band_get_min_value(), rt_band_get_nodata(), rt_band_get_pixtype(), rt_pixtype_index_from_name(), rt_pixtype_name(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_is_empty(), rt_raster_iterator(), rt_raster_new(), rt_raster_serialize(), rt_util_extent_type(), rtpg_nmapalgebra_arg_destroy(), rtpg_nmapalgebra_arg_init(), rtpg_nmapalgebra_callback(), rtpg_nmapalgebra_rastbandarg_process(), rtpg_strtoupper(), rtpg_trim(), rt_raster_serialized_t::size, TRUE, rtpg_nmapalgebra_callback_arg::ufc_info, rtpg_nmapalgebra_callback_arg::ufc_noid, rtpg_nmapalgebra_callback_arg::ufc_rettype, rtpg_nmapalgebra_callback_arg::ufl_info, rt_mask_t::values, rt_mask_t::weighted, pixval::x, and pixval::y.

Referenced by rtpg_nmapalgebra_callback().

509 {
510  rtpg_nmapalgebra_arg arg = NULL;
511  rt_iterator itrset;
512  ArrayType *maskArray;
513  Oid etype;
514  Datum *maskElements;
515  bool *maskNulls;
516  int16 typlen;
517  bool typbyval;
518  char typalign;
519  int ndims = 0;
520  int num;
521  int *maskDims;
522  int x,y;
523 
524 
525  int i = 0;
526  int noerr = 0;
527  int allnull = 0;
528  int allempty = 0;
529  int noband = 0;
530 
531  rt_raster raster = NULL;
532  rt_band band = NULL;
533  rt_pgraster *pgraster = NULL;
534 
535  POSTGIS_RT_DEBUG(3, "Starting...");
536 
537  if (PG_ARGISNULL(0))
538  PG_RETURN_NULL();
539 
540  /* init argument struct */
542  if (arg == NULL) {
543  elog(ERROR, "RASTER_nMapAlgebra: Could not initialize argument structure");
544  PG_RETURN_NULL();
545  }
546 
547  /* let helper function process rastbandarg (0) */
548  if (!rtpg_nmapalgebra_rastbandarg_process(arg, PG_GETARG_ARRAYTYPE_P(0), &allnull, &allempty, &noband)) {
550  elog(ERROR, "RASTER_nMapAlgebra: Could not process rastbandarg");
551  PG_RETURN_NULL();
552  }
553 
554  POSTGIS_RT_DEBUGF(4, "allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
555 
556  /* all rasters are NULL, return NULL */
557  if (allnull == arg->numraster) {
558  elog(NOTICE, "All input rasters are NULL. Returning NULL");
560  PG_RETURN_NULL();
561  }
562 
563  /* pixel type (2) */
564  if (!PG_ARGISNULL(2)) {
565  char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
566 
567  /* Get the pixel type index */
568  arg->pixtype = rt_pixtype_index_from_name(pixtypename);
569  if (arg->pixtype == PT_END) {
571  elog(ERROR, "RASTER_nMapAlgebra: Invalid pixel type: %s", pixtypename);
572  PG_RETURN_NULL();
573  }
574  }
575 
576  /* distancex (3) */
577  if (!PG_ARGISNULL(3)){
578  arg->distance[0] = PG_GETARG_INT32(3);
579  }else{
580  arg->distance[0] = 0;
581  }
582  /* distancey (4) */
583  if (!PG_ARGISNULL(4)){
584  arg->distance[1] = PG_GETARG_INT32(4);
585  }else{
586  arg->distance[1] = 0;
587  }
588  if (arg->distance[0] < 0 || arg->distance[1] < 0) {
590  elog(ERROR, "RASTER_nMapAlgebra: Distance for X and Y axis must be greater than or equal to zero");
591  PG_RETURN_NULL();
592  }
593 
594  /* extent type (5) */
595  if (!PG_ARGISNULL(5)) {
596  char *extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(5))));
597  arg->extenttype = rt_util_extent_type(extenttypename);
598  }
599  POSTGIS_RT_DEBUGF(4, "extenttype: %d", arg->extenttype);
600 
601  /* custom extent (6) */
602  if (arg->extenttype == ET_CUSTOM) {
603  if (PG_ARGISNULL(6)) {
604  elog(NOTICE, "Custom extent is NULL. Returning NULL");
606  PG_RETURN_NULL();
607  }
608 
609  arg->pgcextent = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(6));
610 
611  /* only need the raster header */
613  if (arg->cextent == NULL) {
615  elog(ERROR, "RASTER_nMapAlgebra: Could not deserialize custom extent");
616  PG_RETURN_NULL();
617  }
618  else if (rt_raster_is_empty(arg->cextent)) {
619  elog(NOTICE, "Custom extent is an empty raster. Returning empty raster");
621 
622  raster = rt_raster_new(0, 0);
623  if (raster == NULL) {
624  elog(ERROR, "RASTER_nMapAlgebra: Could not create empty raster");
625  PG_RETURN_NULL();
626  }
627 
628  pgraster = rt_raster_serialize(raster);
629  rt_raster_destroy(raster);
630  if (!pgraster) PG_RETURN_NULL();
631 
632  SET_VARSIZE(pgraster, pgraster->size);
633  PG_RETURN_POINTER(pgraster);
634  }
635  }
636 
637  /* mask (7) */
638 
639  if( PG_ARGISNULL(7) ){
640  pfree(arg->mask);
641  arg->mask = NULL;
642  }
643  else {
644  maskArray = PG_GETARG_ARRAYTYPE_P(7);
645  etype = ARR_ELEMTYPE(maskArray);
646  get_typlenbyvalalign(etype,&typlen,&typbyval,&typalign);
647 
648  switch (etype) {
649  case FLOAT4OID:
650  case FLOAT8OID:
651  break;
652  default:
654  elog(ERROR,"RASTER_nMapAlgebra: Mask data type must be FLOAT8 or FLOAT4");
655  PG_RETURN_NULL();
656  }
657 
658  ndims = ARR_NDIM(maskArray);
659 
660  if (ndims != 2) {
661  elog(ERROR, "RASTER_nMapAlgebra: Mask Must be a 2D array");
663  PG_RETURN_NULL();
664  }
665 
666  maskDims = ARR_DIMS(maskArray);
667 
668  if (maskDims[0] % 2 == 0 || maskDims[1] % 2 == 0) {
669  elog(ERROR,"RASTER_nMapAlgebra: Mask dimensions must be odd");
671  PG_RETURN_NULL();
672  }
673 
674  deconstruct_array(
675  maskArray,
676  etype,
677  typlen, typbyval,typalign,
678  &maskElements,&maskNulls,&num
679  );
680 
681  if (num < 1 || num != (maskDims[0] * maskDims[1])) {
682  if (num) {
683  pfree(maskElements);
684  pfree(maskNulls);
685  }
686  elog(ERROR, "RASTER_nMapAlgebra: Could not deconstruct new values array");
688  PG_RETURN_NULL();
689  }
690 
691  /* allocate mem for mask array */
692  arg->mask->values = palloc(sizeof(double*)* maskDims[0]);
693  arg->mask->nodata = palloc(sizeof(int*)*maskDims[0]);
694  for (i = 0; i < maskDims[0]; i++) {
695  arg->mask->values[i] = (double*) palloc(sizeof(double) * maskDims[1]);
696  arg->mask->nodata[i] = (int*) palloc(sizeof(int) * maskDims[1]);
697  }
698 
699  /* place values in to mask */
700  i = 0;
701  for (y = 0; y < maskDims[0]; y++) {
702  for (x = 0; x < maskDims[1]; x++) {
703  if (maskNulls[i]) {
704  arg->mask->values[y][x] = 0;
705  arg->mask->nodata[y][x] = 1;
706  }
707  else {
708  switch (etype) {
709  case FLOAT4OID:
710  arg->mask->values[y][x] = (double) DatumGetFloat4(maskElements[i]);
711  arg->mask->nodata[y][x] = 0;
712  break;
713  case FLOAT8OID:
714  arg->mask->values[y][x] = (double) DatumGetFloat8(maskElements[i]);
715  arg->mask->nodata[y][x] = 0;
716  }
717  }
718  i++;
719  }
720  }
721 
722  /*set mask dimensions*/
723  arg->mask->dimx = maskDims[0];
724  arg->mask->dimy = maskDims[1];
725  if (maskDims[0] == 1 && maskDims[1] == 1) {
726  arg->distance[0] = 0;
727  arg->distance[1] = 0;
728  }
729  else {
730  arg->distance[0] = maskDims[0] % 2;
731  arg->distance[1] = maskDims[1] % 2;
732  }
733  }/*end if else argisnull*/
734 
735  /* (8) weighted boolean */
736  if (PG_ARGISNULL(8) || !PG_GETARG_BOOL(8)) {
737  if (arg->mask != NULL)
738  arg->mask->weighted = 0;
739  }else{
740  if(arg->mask !=NULL)
741  arg->mask->weighted = 1;
742  }
743 
744  noerr = 1;
745 
746  /* all rasters are empty, return empty raster */
747  if (allempty == arg->numraster) {
748  elog(NOTICE, "All input rasters are empty. Returning empty raster");
749  noerr = 0;
750  }
751  /* all rasters don't have indicated band, return empty raster */
752  else if (noband == arg->numraster) {
753  elog(NOTICE, "All input rasters do not have bands at indicated indexes. Returning empty raster");
754  noerr = 0;
755  }
756  if (!noerr) {
758 
759  raster = rt_raster_new(0, 0);
760  if (raster == NULL) {
761  elog(ERROR, "RASTER_nMapAlgebra: Could not create empty raster");
762  PG_RETURN_NULL();
763  }
764 
765  pgraster = rt_raster_serialize(raster);
766  rt_raster_destroy(raster);
767  if (!pgraster) PG_RETURN_NULL();
768 
769  SET_VARSIZE(pgraster, pgraster->size);
770  PG_RETURN_POINTER(pgraster);
771  }
772 
773  /* do regprocedure last (1) */
774  if (!PG_ARGISNULL(1) || get_fn_expr_argtype(fcinfo->flinfo, 1) == REGPROCEDUREOID) {
775  POSTGIS_RT_DEBUG(4, "processing callbackfunc");
776  arg->callback.ufc_noid = PG_GETARG_OID(1);
777 
778  /* get function info */
779  fmgr_info(arg->callback.ufc_noid, &(arg->callback.ufl_info));
780 
781  /* function cannot return set */
782  noerr = 0;
783  if (arg->callback.ufl_info.fn_retset) {
784  noerr = 1;
785  }
786  /* function should have correct # of args */
787  else if (arg->callback.ufl_info.fn_nargs != 3) {
788  noerr = 2;
789  }
790 
791  /* check that callback function return type is supported */
792  if (
793  get_func_result_type(
794  arg->callback.ufc_noid,
795  &(arg->callback.ufc_rettype),
796  NULL
797  ) != TYPEFUNC_SCALAR
798  ) {
799  noerr = 3;
800  }
801 
802  if (!(
803  arg->callback.ufc_rettype == FLOAT8OID ||
804  arg->callback.ufc_rettype == FLOAT4OID ||
805  arg->callback.ufc_rettype == INT4OID ||
806  arg->callback.ufc_rettype == INT2OID
807  )) {
808  noerr = 4;
809  }
810 
811  /*
812  TODO: consider adding checks of the userfunction parameters
813  should be able to use get_fn_expr_argtype() of fmgr.c
814  */
815 
816  if (noerr != 0) {
818  switch (noerr) {
819  case 4:
820  elog(ERROR, "RASTER_nMapAlgebra: Function provided must return a double precision, float, int or smallint");
821  break;
822  case 3:
823  elog(ERROR, "RASTER_nMapAlgebra: Function provided must return scalar (double precision, float, int, smallint)");
824  break;
825  case 2:
826  elog(ERROR, "RASTER_nMapAlgebra: Function provided must have three input parameters");
827  break;
828  case 1:
829  elog(ERROR, "RASTER_nMapAlgebra: Function provided must return double precision, not resultset");
830  break;
831  }
832  PG_RETURN_NULL();
833  }
834 
835  if (func_volatile(arg->callback.ufc_noid) == 'v')
836  elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
837 
838  /* prep function call data */
839 #if POSTGIS_PGSQL_VERSION > 90
840  InitFunctionCallInfoData(arg->callback.ufc_info, &(arg->callback.ufl_info), arg->callback.ufl_info.fn_nargs, InvalidOid, NULL, NULL);
841 #else
842  InitFunctionCallInfoData(arg->callback.ufc_info, &(arg->callback.ufl_info), arg->callback.ufl_info.fn_nargs, NULL, NULL);
843 #endif
844  memset(arg->callback.ufc_info.argnull, FALSE, sizeof(bool) * arg->callback.ufl_info.fn_nargs);
845 
846  /* userargs (7) */
847  if (!PG_ARGISNULL(9))
848  arg->callback.ufc_info.arg[2] = PG_GETARG_DATUM(9);
849  else {
850  if (arg->callback.ufl_info.fn_strict) {
851  /* build and assign an empty TEXT array */
852  /* TODO: manually free the empty array? */
853  arg->callback.ufc_info.arg[2] = PointerGetDatum(
854  construct_empty_array(TEXTOID)
855  );
856  arg->callback.ufc_info.argnull[2] = FALSE;
857  }
858  else {
859  arg->callback.ufc_info.arg[2] = (Datum) NULL;
860  arg->callback.ufc_info.argnull[2] = TRUE;
861  }
862  }
863  }
864  else {
866  elog(ERROR, "RASTER_nMapAlgebra: callbackfunc must be provided");
867  PG_RETURN_NULL();
868  }
869 
870  /* determine nodataval and possibly pixtype */
871  /* band to check */
872  switch (arg->extenttype) {
873  case ET_LAST:
874  i = arg->numraster - 1;
875  break;
876  case ET_SECOND:
877  if (arg->numraster > 1) {
878  i = 1;
879  break;
880  }
881  default:
882  i = 0;
883  break;
884  }
885  /* find first viable band */
886  if (!arg->hasband[i]) {
887  for (i = 0; i < arg->numraster; i++) {
888  if (arg->hasband[i])
889  break;
890  }
891  if (i >= arg->numraster)
892  i = arg->numraster - 1;
893  }
894  band = rt_raster_get_band(arg->raster[i], arg->nband[i]);
895 
896  /* set pixel type if PT_END */
897  if (arg->pixtype == PT_END)
898  arg->pixtype = rt_band_get_pixtype(band);
899 
900  /* set hasnodata and nodataval */
901  arg->hasnodata = 1;
902  if (rt_band_get_hasnodata_flag(band))
903  rt_band_get_nodata(band, &(arg->nodataval));
904  else
905  arg->nodataval = rt_band_get_min_value(band);
906 
907  POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->pixtype), arg->hasnodata, arg->nodataval);
908 
909  /* init itrset */
910  itrset = palloc(sizeof(struct rt_iterator_t) * arg->numraster);
911  if (itrset == NULL) {
913  elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
914  PG_RETURN_NULL();
915  }
916 
917  /* set itrset */
918  for (i = 0; i < arg->numraster; i++) {
919  itrset[i].raster = arg->raster[i];
920  itrset[i].nband = arg->nband[i];
921  itrset[i].nbnodata = 1;
922  }
923 
924  /* pass everything to iterator */
925  noerr = rt_raster_iterator(
926  itrset, arg->numraster,
927  arg->extenttype, arg->cextent,
928  arg->pixtype,
929  arg->hasnodata, arg->nodataval,
930  arg->distance[0], arg->distance[1],
931  arg->mask,
932  &(arg->callback),
934  &raster
935  );
936 
937  /* cleanup */
938  pfree(itrset);
940 
941  if (noerr != ES_NONE) {
942  elog(ERROR, "RASTER_nMapAlgebra: Could not run raster iterator function");
943  PG_RETURN_NULL();
944  }
945  else if (raster == NULL)
946  PG_RETURN_NULL();
947 
948  pgraster = rt_raster_serialize(raster);
949  rt_raster_destroy(raster);
950 
951  POSTGIS_RT_DEBUG(3, "Finished");
952 
953  if (!pgraster)
954  PG_RETURN_NULL();
955 
956  SET_VARSIZE(pgraster, pgraster->size);
957  PG_RETURN_POINTER(pgraster);
958 }
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
Definition: rt_serialize.c:521
static int rtpg_nmapalgebra_callback(rt_iterator_arg arg, void *userarg, double *value, int *nodata)
raster
Be careful!! Zeros function&#39;s input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
band
Definition: ovdump.py:57
int ** nodata
Definition: librtcore.h:2310
int weighted
Definition: librtcore.h:2311
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition: rt_raster.c:1338
uint16_t dimy
Definition: librtcore.h:2308
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:57
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1597
static void rtpg_nmapalgebra_arg_destroy(rtpg_nmapalgebra_arg arg)
double ** values
Definition: librtcore.h:2309
rtpg_nmapalgebra_callback_arg callback
rt_extenttype rt_util_extent_type(const char *name)
Definition: rt_util.c:191
rt_pixtype rt_pixtype_index_from_name(const char *pixname)
Definition: rt_pixel.c:80
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
static int rtpg_nmapalgebra_rastbandarg_process(rtpg_nmapalgebra_arg arg, ArrayType *array, int *allnull, int *allempty, int *noband)
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:541
char * rtpg_trim(const char *input)
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
rt_raster raster
Definition: librtcore.h:2406
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
uint16_t nband
Definition: librtcore.h:2407
#define FALSE
Definition: dbfopen.c:168
uint8_t nbnodata
Definition: librtcore.h:2408
Struct definitions.
Definition: librtcore.h:2213
double rt_band_get_min_value(rt_band band)
Returns the minimal possible value for the band according to the pixel type.
Definition: rt_band.c:1612
const char * rt_pixtype_name(rt_pixtype pixtype)
Definition: rt_pixel.c:110
rt_errorstate rt_raster_iterator(rt_iterator itrset, uint16_t itrcount, rt_extenttype extenttype, rt_raster customextent, rt_pixtype pixtype, uint8_t hasnodata, double nodataval, uint16_t distancex, uint16_t distancey, rt_mask mask, void *userarg, int(*callback)(rt_iterator_arg arg, void *userarg, double *value, int *nodata), rt_raster *rtnraster)
n-raster iterator.
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:498
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rtpostgis.h:53
static rtpg_nmapalgebra_arg rtpg_nmapalgebra_arg_init()
uint16_t dimx
Definition: librtcore.h:2307
#define TRUE
Definition: dbfopen.c:169
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:717
char * rtpg_strtoupper(char *str)
FunctionCallInfoData ufc_info
Here is the call graph for this function:
Here is the caller graph for this function: