PostGIS 3.0.6dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches

◆ 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);
664 rt_raster_destroy(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);
801 rt_raster_destroy(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 i = (arg->numraster > 1) ? 1 : 0;
937 break;
938 default:
939 i = 0;
940 break;
941 }
942 /* find first viable band */
943 if (!arg->hasband[i]) {
944 for (i = 0; i < arg->numraster; i++) {
945 if (arg->hasband[i])
946 break;
947 }
948 if (i >= arg->numraster)
949 i = arg->numraster - 1;
950 }
951 band = rt_raster_get_band(arg->raster[i], arg->nband[i]);
952
953 /* set pixel type if PT_END */
954 if (arg->pixtype == PT_END)
955 arg->pixtype = rt_band_get_pixtype(band);
956
957 /* set hasnodata and nodataval */
958 arg->hasnodata = 1;
960 rt_band_get_nodata(band, &(arg->nodataval));
961 else
962 arg->nodataval = rt_band_get_min_value(band);
963
964 POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->pixtype), arg->hasnodata, arg->nodataval);
965
966 /* init itrset */
967 itrset = palloc(sizeof(struct rt_iterator_t) * arg->numraster);
968 if (itrset == NULL) {
970 elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
971 PG_RETURN_NULL();
972 }
973
974 /* set itrset */
975 for (i = 0; i < arg->numraster; i++) {
976 itrset[i].raster = arg->raster[i];
977 itrset[i].nband = arg->nband[i];
978 itrset[i].nbnodata = 1;
979 }
980
981 /* pass everything to iterator */
982 noerr = rt_raster_iterator(
983 itrset, arg->numraster,
984 arg->extenttype, arg->cextent,
985 arg->pixtype,
986 arg->hasnodata, arg->nodataval,
987 arg->distance[0], arg->distance[1],
988 arg->mask,
989 &(arg->callback),
991 &raster
992 );
993
994 /* cleanup */
995 pfree(itrset);
997
998 if (noerr != ES_NONE) {
999 elog(ERROR, "RASTER_nMapAlgebra: Could not run raster iterator function");
1000 PG_RETURN_NULL();
1001 }
1002 else if (raster == NULL)
1003 PG_RETURN_NULL();
1004
1005 pgraster = rt_raster_serialize(raster);
1006 rt_raster_destroy(raster);
1007
1008 POSTGIS_RT_DEBUG(3, "Finished");
1009
1010 if (!pgraster)
1011 PG_RETURN_NULL();
1012
1013 SET_VARSIZE(pgraster, pgraster->size);
1014 PG_RETURN_POINTER(pgraster);
1015}
#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:193
const char * rt_pixtype_name(rt_pixtype pixtype)
Definition rt_pixel.c:110
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
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition rt_raster.c:1329
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition rt_raster.c:381
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_trim(const char *input)
char * rtpg_strtoupper(char *str)
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:2444
uint16_t nband
Definition librtcore.h:2445
uint8_t nbnodata
Definition librtcore.h:2446
double ** values
Definition librtcore.h:2347
uint16_t dimy
Definition librtcore.h:2346
uint16_t dimx
Definition librtcore.h:2345
int ** nodata
Definition librtcore.h:2348
Struct definitions.
Definition librtcore.h:2251
rtpg_nmapalgebra_callback_arg callback
FunctionCallInfoData ufc_info

References 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, 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, and rt_mask_t::weighted.

Here is the call graph for this function: