PostGIS  3.4.0dev-r@@SVN_REVISION@@

◆ RASTER_InterpolateRaster()

Datum RASTER_InterpolateRaster ( PG_FUNCTION_ARGS  )

Definition at line 642 of file rtpg_gdal.c.

643 {
644  rt_pgraster *in_pgrast = NULL;
645  rt_pgraster *out_pgrast = NULL;
646  rt_raster in_rast = NULL;
647  rt_raster out_rast = NULL;
648  uint32_t out_rast_bands[1] = {0};
649  rt_band in_band = NULL;
650  rt_band out_band = NULL;
651  int band_number;
652  uint16_t in_band_width, in_band_height;
653  uint32_t npoints;
654  rt_pixtype in_band_pixtype;
655  GDALDataType in_band_gdaltype;
656  size_t in_band_gdaltype_size;
657 
658  rt_envelope env;
659 
660  GDALGridAlgorithm algorithm;
661  text *options_txt = NULL;
662  void *options_struct = NULL;
663  CPLErr err;
664  uint8_t *out_data;
665  rt_errorstate rterr;
666 
667  /* Input points */
668  LWPOINTITERATOR *iterator;
669  POINT4D pt;
670  size_t coord_count = 0;
671  LWGEOM *lwgeom;
672  double *xcoords, *ycoords, *zcoords;
673 
674  GSERIALIZED *gser = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
675 
676  /* Z value is required to drive the grid heights */
677  if (!gserialized_has_z(gser))
678  elog(ERROR, "%s: input geometry does not have Z values", __func__);
679 
680  /* Cannot process empties */
681  if (gserialized_is_empty(gser))
682  PG_RETURN_NULL();
683 
684  in_pgrast = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(2));
685  in_rast = rt_raster_deserialize(in_pgrast, FALSE);
686  if (!in_rast)
687  elog(ERROR, "%s: Could not deserialize raster", __func__);
688 
689  /* GDAL cannot grid a skewed raster */
690  if (rt_raster_get_x_skew(in_rast) != 0.0 ||
691  rt_raster_get_y_skew(in_rast) != 0.0) {
692  elog(ERROR, "%s: Cannot generate a grid into a skewed raster",__func__);
693  }
694 
695  /* Flat JSON map of options from user */
696  options_txt = PG_GETARG_TEXT_P(1);
697  /* 1-base band number from user */
698  band_number = PG_GETARG_INT32(3);
699  if (band_number < 1)
700  elog(ERROR, "%s: Invalid band number %d", __func__, band_number);
701 
702  lwgeom = lwgeom_from_gserialized(gser);
703  npoints = lwgeom_count_vertices(lwgeom);
704  /* This shouldn't happen, but just in case... */
705  if (npoints < 1)
706  elog(ERROR, "%s: Geometry has no points", __func__);
707 
708  in_band = rt_raster_get_band(in_rast, band_number-1);
709  if (!in_band)
710  elog(ERROR, "%s: Cannot access raster band %d", __func__, band_number);
711 
712 
713  rterr = rt_raster_get_envelope(in_rast, &env);
714  if (rterr == ES_ERROR)
715  elog(ERROR, "%s: Unable to calculate envelope", __func__);
716 
717  /* Get geometry of input raster */
718  in_band_width = rt_band_get_width(in_band);
719  in_band_height = rt_band_get_height(in_band);
720  in_band_pixtype = rt_band_get_pixtype(in_band);
721  in_band_gdaltype = rt_util_pixtype_to_gdal_datatype(in_band_pixtype);
722  in_band_gdaltype_size = GDALGetDataTypeSize(in_band_gdaltype) / 8;
723 
724  /* Quickly copy options struct into local memory context, so we */
725  /* don't have malloc'ed memory lying around */
726  // if (err == CE_None && options_struct) {
727  // void *tmp = options_struct;
728  // switch (algorithm) {
729  // case GGA_InverseDistanceToAPower:
730  // options_struct = palloc(sizeof(GDALGridInverseDistanceToAPowerOptions));
731  // memcpy(options_struct, tmp, sizeof(GDALGridInverseDistanceToAPowerOptions));
732  // break;
733  // case GGA_InverseDistanceToAPowerNearestNeighbor:
734  // options_struct = palloc(sizeof(GDALGridInverseDistanceToAPowerNearestNeighborOptions));
735  // memcpy(options_struct, tmp, sizeof(GDALGridInverseDistanceToAPowerNearestNeighborOptions));
736  // break;
737  // case GGA_MovingAverage:
738  // options_struct = palloc(sizeof(GDALGridMovingAverageOptions));
739  // memcpy(options_struct, tmp, sizeof(GDALGridMovingAverageOptions));
740  // break;
741  // case GGA_NearestNeighbor:
742  // options_struct = palloc(sizeof(GDALGridNearestNeighborOptions));
743  // memcpy(options_struct, tmp, sizeof(GDALGridNearestNeighborOptions));
744  // break;
745  // case GGA_Linear:
746  // options_struct = palloc(sizeof(GDALGridLinearOptions));
747  // memcpy(options_struct, tmp, sizeof(GDALGridLinearOptions));
748  // break;
749  // default:
750  // elog(ERROR, "%s: Unsupported gridding algorithm %d", __func__, algorithm);
751  // }
752  // free(tmp);
753  // }
754 
755  /* Prepare destination grid buffer for output */
756  out_data = palloc(in_band_gdaltype_size * in_band_width * in_band_height);
757 
758  /* Prepare input points for processing */
759  xcoords = palloc(sizeof(double) * npoints);
760  ycoords = palloc(sizeof(double) * npoints);
761  zcoords = palloc(sizeof(double) * npoints);
762 
763  /* Populate input points */
764  iterator = lwpointiterator_create(lwgeom);
765  while(lwpointiterator_next(iterator, &pt) == LW_SUCCESS) {
766  if (coord_count >= npoints)
767  elog(ERROR, "%s: More points from iterator than expected", __func__);
768  xcoords[coord_count] = pt.x;
769  ycoords[coord_count] = pt.y;
770  zcoords[coord_count] = pt.z;
771  coord_count++;
772  }
773  lwpointiterator_destroy(iterator);
774 
775  /* Extract algorithm and options from options text */
776  /* This malloc's the options struct, so clean up right away */
777  err = ParseAlgorithmAndOptions(
778  text_to_cstring(options_txt),
779  &algorithm,
780  &options_struct);
781  if (err != CE_None) {
782  if (options_struct) free(options_struct);
783  elog(ERROR, "%s: Unable to parse options string: %s", __func__, CPLGetLastErrorMsg());
784  }
785 
786  /* Run the gridding algorithm */
787  err = GDALGridCreate(
788  algorithm, options_struct,
789  npoints, xcoords, ycoords, zcoords,
790  env.MinX, env.MaxX, env.MinY, env.MaxY,
791  in_band_width, in_band_height,
792  in_band_gdaltype, out_data,
793  NULL, /* GDALProgressFunc */
794  NULL /* ProgressArgs */
795  );
796 
797  /* Quickly clean up malloc'ed memory */
798  if (options_struct)
799  free(options_struct);
800 
801  if (err != CE_None) {
802  elog(ERROR, "%s: GDALGridCreate failed: %s", __func__, CPLGetLastErrorMsg());
803  }
804 
805  out_rast_bands[0] = band_number-1;
806  out_rast = rt_raster_from_band(in_rast, out_rast_bands, 1);
807  out_band = rt_raster_get_band(out_rast, 0);
808  if (!out_band)
809  elog(ERROR, "%s: Cannot access output raster band", __func__);
810 
811  /* Copy the data from the output buffer into the destination band */
812  for (uint32_t y = 0; y < in_band_height; y++) {
813  size_t offset = (in_band_height-y-1) * (in_band_gdaltype_size * in_band_width);
814  rterr = rt_band_set_pixel_line(out_band, 0, y, out_data + offset, in_band_width);
815  }
816 
817  out_pgrast = rt_raster_serialize(out_rast);
818  rt_raster_destroy(out_rast);
819  rt_raster_destroy(in_rast);
820 
821  if (NULL == out_pgrast) PG_RETURN_NULL();
822 
823  SET_VARSIZE(out_pgrast, out_pgrast->size);
824  PG_RETURN_POINTER(out_pgrast);
825 }
#define FALSE
Definition: dbfopen.c:72
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
Definition: gserialized.c:239
int gserialized_is_empty(const GSERIALIZED *g)
Check if a GSERIALIZED is empty without deserializing first.
Definition: gserialized.c:152
int gserialized_has_z(const GSERIALIZED *g)
Check if a GSERIALIZED has a Z ordinate.
Definition: gserialized.c:174
LWPOINTITERATOR * lwpointiterator_create(const LWGEOM *g)
Create a new LWPOINTITERATOR over supplied LWGEOM*.
Definition: lwiterator.c:242
int lwpointiterator_next(LWPOINTITERATOR *s, POINT4D *p)
Attempts to assign the next point in the iterator to p, and advances the iterator to the next point.
Definition: lwiterator.c:210
#define LW_SUCCESS
Definition: liblwgeom.h:97
void lwpointiterator_destroy(LWPOINTITERATOR *s)
Free all memory associated with the iterator.
Definition: lwiterator.c:267
uint32_t lwgeom_count_vertices(const LWGEOM *geom)
Count the total number of vertices in any LWGEOM.
Definition: lwgeom.c:1246
uint16_t rt_band_get_width(rt_band band)
Return width of this band.
Definition: rt_band.c:640
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
Definition: rt_raster.c:185
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:86
rt_pixtype
Definition: librtcore.h:187
GDALDataType rt_util_pixtype_to_gdal_datatype(rt_pixtype pt)
Convert rt_pixtype to GDALDataType.
Definition: rt_util.c:124
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
Definition: rt_serialize.c:521
rt_errorstate
Enum definitions.
Definition: librtcore.h:181
@ ES_ERROR
Definition: librtcore.h:183
rt_errorstate rt_band_set_pixel_line(rt_band band, int x, int y, void *vals, uint32_t len)
Set values of multiple pixels.
Definition: rt_band.c:853
rt_raster rt_raster_from_band(rt_raster raster, uint32_t *bandNums, int count)
Construct a new rt_raster from an existing rt_raster and an array of band numbers.
Definition: rt_raster.c:1463
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:631
rt_errorstate rt_raster_get_envelope(rt_raster raster, rt_envelope *env)
Get raster's envelope.
Definition: rt_raster.c:904
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
Definition: rt_raster.c:194
uint16_t rt_band_get_height(rt_band band)
Return height of this band.
Definition: rt_band.c:649
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:725
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:385
void free(void *)
double x
Definition: liblwgeom.h:414
double z
Definition: liblwgeom.h:414
double y
Definition: liblwgeom.h:414
double MinX
Definition: librtcore.h:167
double MaxX
Definition: librtcore.h:168
double MinY
Definition: librtcore.h:169
double MaxY
Definition: librtcore.h:170
Struct definitions.
Definition: librtcore.h:2403

References ES_ERROR, FALSE, free(), gserialized_has_z(), gserialized_is_empty(), LW_SUCCESS, lwgeom_count_vertices(), lwgeom_from_gserialized(), lwpointiterator_create(), lwpointiterator_destroy(), lwpointiterator_next(), rt_envelope::MaxX, rt_envelope::MaxY, rt_envelope::MinX, rt_envelope::MinY, rt_band_get_height(), rt_band_get_pixtype(), rt_band_get_width(), rt_band_set_pixel_line(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_from_band(), rt_raster_get_band(), rt_raster_get_envelope(), rt_raster_get_x_skew(), rt_raster_get_y_skew(), rt_raster_serialize(), rt_util_pixtype_to_gdal_datatype(), rt_raster_serialized_t::size, POINT4D::x, POINT4D::y, pixval::y, and POINT4D::z.

Here is the call graph for this function: