PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ RASTER_nMapAlgebra()

Datum RASTER_nMapAlgebra ( PG_FUNCTION_ARGS  )

Definition at line 543 of file rtpg_mapalgebra.c.

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

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, rt_iterator_t::nband, rtpg_nmapalgebra_arg_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, rt_iterator_t::raster, rtpg_nmapalgebra_arg_t::raster, rtrowdump::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, text_to_cstring(), 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.

Here is the call graph for this function: