RASTER_InterpolateRaster()

Datum RASTER_InterpolateRaster

(

PG_FUNCTION_ARGS

)

Definition at line 643 of file rtpg_gdal.c.

{
        rt_pgraster *in_pgrast = NULL;
        rt_pgraster *out_pgrast = NULL;
        rt_raster in_rast = NULL;
        rt_raster out_rast = NULL;
        uint32_t out_rast_bands[1] = {0};
        rt_band in_band = NULL;
        rt_band out_band = NULL;
        int band_number;
        uint16_t in_band_width, in_band_height;
        uint32_t npoints;
        rt_pixtype in_band_pixtype;
        GDALDataType in_band_gdaltype;
        size_t in_band_gdaltype_size;
 
        rt_envelope env;
 
        GDALGridAlgorithm algorithm;
        text *options_txt = NULL;
        void *options_struct = NULL;
        CPLErr err;
        uint8_t *out_data;
        rt_errorstate rterr;
 
        /* Input points */
        LWPOINTITERATOR *iterator;
        POINT4D pt;
        size_t coord_count = 0;
        LWGEOM *lwgeom;
        double *xcoords, *ycoords, *zcoords;
 
        GSERIALIZED *gser = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
 
        /* Z value is required to drive the grid heights */
        if (!gserialized_has_z(gser))
                elog(ERROR, "%s: input geometry does not have Z values", __func__);
 
        /* Cannot process empties */
        if (gserialized_is_empty(gser))
                PG_RETURN_NULL();
 
        in_pgrast = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(2));
        in_rast = rt_raster_deserialize(in_pgrast, FALSE);
        if (!in_rast)
                elog(ERROR, "%s: Could not deserialize raster", __func__);
 
        /* GDAL cannot grid a skewed raster */
        if (rt_raster_get_x_skew(in_rast) != 0.0 ||
            rt_raster_get_y_skew(in_rast) != 0.0) {
                elog(ERROR, "%s: Cannot generate a grid into a skewed raster",__func__);
        }
 
        /* Flat JSON map of options from user */
        options_txt = PG_GETARG_TEXT_P(1);
        /* 1-base band number from user */
        band_number = PG_GETARG_INT32(3);
        if (band_number < 1)
                elog(ERROR, "%s: Invalid band number %d", __func__, band_number);
 
        lwgeom = lwgeom_from_gserialized(gser);
        npoints = lwgeom_count_vertices(lwgeom);
        /* This shouldn't happen, but just in case... */
        if (npoints < 1)
                elog(ERROR, "%s: Geometry has no points", __func__);
 
        in_band = rt_raster_get_band(in_rast, band_number-1);
        if (!in_band)
                elog(ERROR, "%s: Cannot access raster band %d", __func__, band_number);
 
 
        rterr = rt_raster_get_envelope(in_rast, &env);
        if (rterr == ES_ERROR)
                elog(ERROR, "%s: Unable to calculate envelope", __func__);
 
        /* Get geometry of input raster */
        in_band_width = rt_band_get_width(in_band);
        in_band_height = rt_band_get_height(in_band);
        in_band_pixtype = rt_band_get_pixtype(in_band);
        in_band_gdaltype = rt_util_pixtype_to_gdal_datatype(in_band_pixtype);
        in_band_gdaltype_size = GDALGetDataTypeSize(in_band_gdaltype) / 8;
 
        /* Quickly copy options struct into local memory context, so we */
        /* don't have malloc'ed memory lying around */
        // if (err == CE_None && options_struct) {
        //      void *tmp = options_struct;
        //      switch (algorithm) {
        //              case GGA_InverseDistanceToAPower:
        //                      options_struct = palloc(sizeof(GDALGridInverseDistanceToAPowerOptions));
        //                      memcpy(options_struct, tmp, sizeof(GDALGridInverseDistanceToAPowerOptions));
        //                      break;
        //              case GGA_InverseDistanceToAPowerNearestNeighbor:
        //                      options_struct = palloc(sizeof(GDALGridInverseDistanceToAPowerNearestNeighborOptions));
        //                      memcpy(options_struct, tmp, sizeof(GDALGridInverseDistanceToAPowerNearestNeighborOptions));
        //                      break;
        //              case GGA_MovingAverage:
        //                      options_struct = palloc(sizeof(GDALGridMovingAverageOptions));
        //                      memcpy(options_struct, tmp, sizeof(GDALGridMovingAverageOptions));
        //                      break;
        //              case GGA_NearestNeighbor:
        //                      options_struct = palloc(sizeof(GDALGridNearestNeighborOptions));
        //                      memcpy(options_struct, tmp, sizeof(GDALGridNearestNeighborOptions));
        //                      break;
        //              case GGA_Linear:
        //                      options_struct = palloc(sizeof(GDALGridLinearOptions));
        //                      memcpy(options_struct, tmp, sizeof(GDALGridLinearOptions));
        //                      break;
        //              default:
        //                      elog(ERROR, "%s: Unsupported gridding algorithm %d", __func__, algorithm);
        //      }
        //      free(tmp);
        // }
 
        /* Prepare destination grid buffer for output */
        out_data = palloc(in_band_gdaltype_size * in_band_width * in_band_height);
 
        /* Prepare input points for processing */
        xcoords = palloc(sizeof(double) * npoints);
        ycoords = palloc(sizeof(double) * npoints);
        zcoords = palloc(sizeof(double) * npoints);
 
        /* Populate input points */
        iterator = lwpointiterator_create(lwgeom);
        while(lwpointiterator_next(iterator, &pt) == LW_SUCCESS) {
                if (coord_count >= npoints)
                        elog(ERROR, "%s: More points from iterator than expected", __func__);
                xcoords[coord_count] = pt.x;
                ycoords[coord_count] = pt.y;
                zcoords[coord_count] = pt.z;
                coord_count++;
        }
        lwpointiterator_destroy(iterator);
 
        /* Extract algorithm and options from options text */
        /* This malloc's the options struct, so clean up right away */
        err = ParseAlgorithmAndOptions(
                text_to_cstring(options_txt),
                &algorithm,
                &options_struct);
        if (err != CE_None) {
                if (options_struct) free(options_struct);
                elog(ERROR, "%s: Unable to parse options string: %s", __func__, CPLGetLastErrorMsg());
        }
 
        /* Run the gridding algorithm */
        err = GDALGridCreate(
                algorithm, options_struct,
                npoints, xcoords, ycoords, zcoords,
                env.MinX, env.MaxX, env.MinY, env.MaxY,
                in_band_width, in_band_height,
                in_band_gdaltype, out_data,
                NULL, /* GDALProgressFunc */
                NULL /* ProgressArgs */
                );
 
        /* Quickly clean up malloc'ed memory */
        if (options_struct)
                free(options_struct);
 
        if (err != CE_None) {
                elog(ERROR, "%s: GDALGridCreate failed: %s", __func__, CPLGetLastErrorMsg());
        }
 
        out_rast_bands[0] = band_number-1;
        out_rast = rt_raster_from_band(in_rast, out_rast_bands, 1);
        out_band = rt_raster_get_band(out_rast, 0);
        if (!out_band)
                elog(ERROR, "%s: Cannot access output raster band", __func__);
 
        /* Copy the data from the output buffer into the destination band */
        for (uint32_t y = 0; y < in_band_height; y++) {
                size_t offset = (in_band_height-y-1) * (in_band_gdaltype_size * in_band_width);
                rterr = rt_band_set_pixel_line(out_band, 0, y, out_data + offset, in_band_width);
        }
 
        out_pgrast = rt_raster_serialize(out_rast);
        rt_raster_destroy(out_rast);
        rt_raster_destroy(in_rast);
 
        if (NULL == out_pgrast) PG_RETURN_NULL();
 
        SET_VARSIZE(out_pgrast, out_pgrast->size);
        PG_RETURN_POINTER(out_pgrast);
}

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: