PostGIS  2.4.9dev-r@@SVN_REVISION@@

◆ RASTER_nMapAlgebra()

Datum RASTER_nMapAlgebra ( PG_FUNCTION_ARGS  )

Definition at line 534 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().

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