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

◆ RASTER_getPixelCentroids()

Datum RASTER_getPixelCentroids ( PG_FUNCTION_ARGS  )

Definition at line 664 of file rtpg_geometry.c.

665{
666 FuncCallContext *funcctx;
667 TupleDesc tupdesc;
668 rt_pixel pix = NULL;
669 rt_pixel pix2;
670 int call_cntr;
671 int max_calls;
672 int i = 0;
673
674 if (SRF_IS_FIRSTCALL()) {
675 MemoryContext oldcontext;
676 rt_pgraster *pgraster = NULL;
677 rt_raster raster = NULL;
678 rt_band band = NULL;
679 int nband = 1;
680 int numbands;
681 bool hasband = FALSE;
682 bool exclude_nodata_value = TRUE;
683 bool nocolumnx = TRUE;
684 bool norowy = TRUE;
685 int x = 0;
686 int y = 0;
687 int bounds[4] = {0};
688 int pixcount = 0;
689 double value = 0;
690 int isnodata = 0;
691
692 LWPOINT* point = NULL;
693
694 POSTGIS_RT_DEBUG(3, "RASTER_getPixelCentroids first call");
695
696 /* create a function context for cross-call persistence */
697 funcctx = SRF_FIRSTCALL_INIT();
698
699 /* switch to memory context appropriate for multiple function calls */
700 oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
701
702 /* raster */
703 if (PG_ARGISNULL(0)) {
704 MemoryContextSwitchTo(oldcontext);
705 SRF_RETURN_DONE(funcctx);
706 }
707 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
708
709 /* band */
710 if (!PG_ARGISNULL(1)) {
711 nband = PG_GETARG_INT32(1);
712 hasband = TRUE;
713 }
714
715 /* column */
716 if (!PG_ARGISNULL(2)) {
717 bounds[0] = PG_GETARG_INT32(2);
718 bounds[1] = bounds[0];
719 nocolumnx = FALSE;
720 }
721
722 /* row */
723 if (!PG_ARGISNULL(3)) {
724 bounds[2] = PG_GETARG_INT32(3);
725 bounds[3] = bounds[2];
726 norowy = FALSE;
727 }
728
729 /* exclude NODATA */
730 if (!PG_ARGISNULL(4))
731 exclude_nodata_value = PG_GETARG_BOOL(4);
732
733 /* deserialize raster */
734 raster = rt_raster_deserialize(pgraster, FALSE);
735 if (!raster) {
736 PG_FREE_IF_COPY(pgraster, 0);
737 ereport(ERROR, (
738 errcode(ERRCODE_OUT_OF_MEMORY),
739 errmsg("Could not deserialize raster")
740 ));
741 MemoryContextSwitchTo(oldcontext);
742 SRF_RETURN_DONE(funcctx);
743 }
744
745 /* raster empty, return empty set */
746 if (rt_raster_is_empty(raster)) {
747 elog(NOTICE, "Raster is empty. Returning empty set");
748 rt_raster_destroy(raster);
749 PG_FREE_IF_COPY(pgraster, 0);
750 MemoryContextSwitchTo(oldcontext);
751 SRF_RETURN_DONE(funcctx);
752 }
753
754 /* band specified, load band and info */
755 if (hasband) {
756 numbands = rt_raster_get_num_bands(raster);
757 POSTGIS_RT_DEBUGF(3, "band %d", nband);
758 POSTGIS_RT_DEBUGF(3, "# of bands %d", numbands);
759
760 /* check band index */
761 if (nband < 1 || nband > numbands) {
762 elog(NOTICE, "Invalid band index (must use 1-based). Returning empty set");
763 rt_raster_destroy(raster);
764 PG_FREE_IF_COPY(pgraster, 0);
765 MemoryContextSwitchTo(oldcontext);
766 SRF_RETURN_DONE(funcctx);
767 }
768
769 /* get band */
770 band = rt_raster_get_band(raster, nband - 1);
771 if (!band) {
772 elog(NOTICE, "Could not find band at index %d. Returning empty set", nband);
773
774 rt_raster_destroy(raster);
775 PG_FREE_IF_COPY(pgraster, 0);
776
777 MemoryContextSwitchTo(oldcontext);
778 SRF_RETURN_DONE(funcctx);
779 }
780
781 if (!rt_band_get_hasnodata_flag(band))
782 exclude_nodata_value = FALSE;
783 }
784
785 /* set bounds if columnx, rowy not set */
786 if (nocolumnx) {
787 bounds[0] = 1;
788 bounds[1] = rt_raster_get_width(raster);
789 }
790 if (norowy) {
791 bounds[2] = 1;
792 bounds[3] = rt_raster_get_height(raster);
793 }
794 POSTGIS_RT_DEBUGF(3, "bounds (min x, max x, min y, max y) = (%d, %d, %d, %d)",
795 bounds[0], bounds[1], bounds[2], bounds[3]);
796
797 /* rowy */
798 pixcount = 0;
799 for (y = bounds[2]; y <= bounds[3]; y++) {
800 /* columnx */
801 for (x = bounds[0]; x <= bounds[1]; x++) {
802
803 value = 0;
804 isnodata = TRUE;
805
806 if (hasband) {
807 if (rt_band_get_pixel(band, x - 1, y - 1, &value, &isnodata) != ES_NONE) {
808 /* free already created centroid points */
809 for (i = 0; i < pixcount; i++)
810 lwgeom_free(pix[i].geom);
811 if (pixcount) pfree(pix);
812
813 rt_band_destroy(band);
814 rt_raster_destroy(raster);
815 PG_FREE_IF_COPY(pgraster, 0);
816
817 MemoryContextSwitchTo(oldcontext);
818 elog(ERROR, "RASTER_getPixelCentroids: Could not get pixel value");
819 SRF_RETURN_DONE(funcctx);
820 }
821
822 /* don't continue if pixel is NODATA and to exclude NODATA */
823 if (isnodata && exclude_nodata_value) {
824 POSTGIS_RT_DEBUG(5, "pixel value is NODATA and exclude_nodata_value = TRUE");
825 continue;
826 }
827 }
828
829 /* geometry */
830 point = rt_raster_pixel_as_centroid_point(raster, x - 1, y - 1);
831 if (!point) {
832 /* free already created centroid points */
833 for (i = 0; i < pixcount; i++)
834 lwgeom_free(pix[i].geom);
835 if (pixcount) pfree(pix);
836
837 if (hasband) rt_band_destroy(band);
838 rt_raster_destroy(raster);
839 PG_FREE_IF_COPY(pgraster, 0);
840
841 MemoryContextSwitchTo(oldcontext);
842 elog(ERROR, "RASTER_getPixelCentroids: Could not get pixel centroid");
843 SRF_RETURN_DONE(funcctx);
844 }
845
846 /* allocate space for new point */
847 if (!pixcount) {
848 pix = palloc(sizeof(struct rt_pixel_t) * (pixcount + 1));
849 }
850 else {
851 pix = repalloc(pix, sizeof(struct rt_pixel_t) * (pixcount + 1));
852 }
853 if (pix == NULL) {
854 lwpoint_free(point);
855
856 if (hasband) rt_band_destroy(band);
857 rt_raster_destroy(raster);
858 PG_FREE_IF_COPY(pgraster, 0);
859
860 MemoryContextSwitchTo(oldcontext);
861 elog(ERROR, "RASTER_getPixelCentroids: Could not allocate memory for storing pixel centroids");
862 SRF_RETURN_DONE(funcctx);
863 }
864
865 /* set pixel geometry */
866 pix[pixcount].geom = (LWGEOM *) point;
867 POSTGIS_RT_DEBUGF(5, "point @ %p", point);
868 POSTGIS_RT_DEBUGF(5, "geom @ %p", pix[pixcount].geom);
869
870 /* set pixel coordinates */
871 pix[pixcount].x = x;
872 pix[pixcount].y = y;
873
874 /* set pixel value */
875 pix[pixcount].value = value;
876
877 /* NODATA */
878 if (hasband) {
879 if (exclude_nodata_value)
880 pix[pixcount].nodata = isnodata;
881 else
882 pix[pixcount].nodata = FALSE;
883 }
884 else {
885 pix[pixcount].nodata = isnodata;
886 }
887
888 /* update pixcount*/
889 pixcount++;
890 }
891 }
892
893 /* cleanup */
894 if (hasband) rt_band_destroy(band);
895 rt_raster_destroy(raster);
896 PG_FREE_IF_COPY(pgraster, 0);
897
898 /* shortcut if no pixcount */
899 if (pixcount < 1) {
900 elog(NOTICE, "No pixels found for band %d", nband);
901 MemoryContextSwitchTo(oldcontext);
902 SRF_RETURN_DONE(funcctx);
903 }
904
905 /* store pixel data */
906 funcctx->user_fctx = pix;
907
908 /* set total number of tuples to be returned */
909 funcctx->max_calls = pixcount;
910 POSTGIS_RT_DEBUGF(3, "pixcount = %d", pixcount);
911
912 /* build a tuple descriptor for our result type */
913 if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
914 ereport(ERROR, (
915 errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
916 errmsg("function returning record called in context that cannot accept type record")
917 ));
918 }
919 /* set tuple descriptor */
920 BlessTupleDesc(tupdesc);
921 funcctx->tuple_desc = tupdesc;
922
923 MemoryContextSwitchTo(oldcontext);
924 }
925
926 /* stuff done on every call of the function */
927 funcctx = SRF_PERCALL_SETUP();
928
929 call_cntr = funcctx->call_cntr;
930 max_calls = funcctx->max_calls;
931 tupdesc = funcctx->tuple_desc;
932 pix2 = funcctx->user_fctx;
933
934 /* do when there is more left to send */
935 if (call_cntr < max_calls) {
936 Datum values[VALUES_LENGTH];
937 bool nulls[VALUES_LENGTH];
938 HeapTuple tuple;
939 Datum result;
940
941 GSERIALIZED *gser = NULL;
942 size_t gser_size = 0;
943
944 POSTGIS_RT_DEBUGF(3, "call number %d", call_cntr);
945
946 memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
947
948 /* convert LWGEOM go GSERIALIZED */
949 gser = gserialized_from_lwgeom(pix2[call_cntr].geom, &gser_size);
950 lwgeom_free(pix2[call_cntr].geom); /* free geometry object */
951
952 /* set result tuple data */
953 values[0] = PointerGetDatum(gser);
954 if (pix2[call_cntr].nodata)
955 nulls[1] = TRUE;
956 else
957 values[1] = Float8GetDatum(pix2[call_cntr].value);
958 values[2] = Int32GetDatum(pix2[call_cntr].x);
959 values[3] = Int32GetDatum(pix2[call_cntr].y);
960
961 /* build a tuple */
962 tuple = heap_form_tuple(tupdesc, values, nulls);
963
964 /* make the tuple into a datum */
965 result = HeapTupleGetDatum(tuple);
966
967 SRF_RETURN_NEXT(funcctx, result);
968 }
969 /* do when there is no more left */
970 else {
971 pfree(pix2);
972 SRF_RETURN_DONE(funcctx);
973 }
974}
result
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition cu_print.c:267
TRUE
#define TRUE
Definition dbfopen.c:73
FALSE
#define FALSE
Definition dbfopen.c:72
gserialized_from_lwgeom
GSERIALIZED * gserialized_from_lwgeom(LWGEOM *geom, size_t *size)
Allocate a new GSERIALIZED from an LWGEOM.
Definition gserialized.c:251
lwpoint_free
void lwpoint_free(LWPOINT *pt)
Definition lwpoint.c:213
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition lwgeom.c:1246
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition rt_band.c:833
rt_band_get_pixel
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition rt_band.c:1551
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition rt_raster.c:86
ES_NONE
@ ES_NONE
Definition librtcore.h:182
rt_band_destroy
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition rt_band.c:499
rt_raster_get_num_bands
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition rt_raster.c:376
rt_raster_get_height
uint16_t rt_raster_get_height(rt_raster raster)
Definition rt_raster.c:133
rt_raster_pixel_as_centroid_point
LWPOINT * rt_raster_pixel_as_centroid_point(rt_raster rast, int x, int y)
Get a raster pixel centroid point.
Definition rt_geometry.c:669
rt_raster_get_width
uint16_t rt_raster_get_width(rt_raster raster)
Definition rt_raster.c:125
rt_raster_deserialize
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition rt_serialize.c:732
rt_raster_is_empty
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition rt_raster.c:1240
rt_raster_get_band
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition rt_raster.c:385
genraster.value
int value
Definition genraster.py:62
pixval.nband
nband
Definition pixval.py:56
rtrowdump.raster
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition rtrowdump.py:125
VALUES_LENGTH
#define VALUES_LENGTH
Definition rtpg_geometry.c:191
POSTGIS_RT_DEBUG
#define POSTGIS_RT_DEBUG(level, msg)
Definition rtpostgis.h:65
POSTGIS_RT_DEBUGF
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition rtpostgis.h:69
rt_band_t
Definition librtcore.h:2515
rt_pixel_t::y
int y
Definition librtcore.h:2537
rt_pixel_t::geom
LWGEOM * geom
Definition librtcore.h:2542
rt_pixel_t::value
double value
Definition librtcore.h:2540
rt_pixel_t::nodata
uint8_t nodata
Definition librtcore.h:2539
rt_pixel_t::x
int x
Definition librtcore.h:2536
rt_pixel_t
Definition librtcore.h:2535
rt_raster_serialized_t
Struct definitions.
Definition librtcore.h:2452
rt_raster_t
Definition librtcore.h:2486

References ES_NONE, FALSE, rt_pixel_t::geom, gserialized_from_lwgeom(), lwgeom_free(), lwpoint_free(), rt_pixel_t::nodata, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, result, rt_band_destroy(), rt_band_get_hasnodata_flag(), rt_band_get_pixel(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_height(), rt_raster_get_num_bands(), rt_raster_get_width(), rt_raster_is_empty(), rt_raster_pixel_as_centroid_point(), TRUE, rt_pixel_t::value, VALUES_LENGTH, rt_pixel_t::x, and rt_pixel_t::y.

Here is the call graph for this function: