PostGIS  2.5.1dev-r@@SVN_REVISION@@

◆ rt_raster_gdal_rasterize()

rt_raster rt_raster_gdal_rasterize ( const unsigned char *  wkb,
uint32_t  wkb_len,
const char *  srs,
uint32_t  num_bands,
rt_pixtype pixtype,
double *  init,
double *  value,
double *  nodata,
uint8_t hasnodata,
int *  width,
int *  height,
double *  scale_x,
double *  scale_y,
double *  ul_xw,
double *  ul_yw,
double *  grid_xw,
double *  grid_yw,
double *  skew_x,
double *  skew_y,
char **  options 
)

Return a raster of the provided geometry.

Parameters
wkb: WKB representation of the geometry to convert
wkb_len: length of the WKB representation of the geometry
srs: the geometry's coordinate system in OGC WKT
num_bands: number of bands in the output raster
pixtype: data type of each band
init: array of values to initialize each band with
value: array of values for pixels of geometry
nodata: array of nodata values for each band
hasnodata: array flagging the presence of nodata for each band
width: the number of columns of the raster
height: the number of rows of the raster
scale_x: the pixel width of the raster
scale_y: the pixel height of the raster
ul_xw: the X value of upper-left corner of the raster
ul_yw: the Y value of upper-left corner of the raster
grid_xw: the X value of point on grid to align raster to
grid_yw: the Y value of point on grid to align raster to
skew_x: the X skew of the raster
skew_y: the Y skew of the raster
options: array of options. only option is "ALL_TOUCHED"
Returns
the raster of the provided geometry or NULL

Definition at line 2509 of file rt_raster.c.

References _rti_rasterize_arg_destroy(), _rti_rasterize_arg_init(), ovdump::band, _rti_rasterize_arg_t::bandlist, ovdump::data, ES_NONE, FLT_EQ, FLT_NEQ, _rti_rasterize_arg_t::hasnodata, rt_raster_t::height, _rti_rasterize_arg_t::init, rt_raster_t::ipX, rt_raster_t::ipY, LW_PARSER_CHECK_NONE, LWGEOM2GEOS(), lwgeom_free(), lwgeom_from_wkb(), lwgeom_geos_error(), lwpoly_as_lwgeom(), lwpoly_free(), rt_envelope::MaxX, rt_envelope::MaxY, rt_envelope::MinX, rt_envelope::MinY, _rti_rasterize_arg_t::noband, _rti_rasterize_arg_t::nodata, _rti_rasterize_arg_t::numbands, _rti_rasterize_arg_t::pixtype, PT_8BUI, rtpixdump::rast, RASTER_DEBUG, RASTER_DEBUGF, rt_band_destroy(), rt_band_get_hasnodata_flag(), rt_band_get_nodata(), rt_band_get_pixel(), rt_band_get_pixtype(), rt_band_new_inline(), rt_band_set_ownsdata_flag(), rt_band_set_pixel(), rt_pixtype_size(), rt_raster_cell_to_geopoint(), rt_raster_compute_skewed_raster(), rt_raster_destroy(), rt_raster_from_gdal_dataset(), rt_raster_geopoint_to_cell(), rt_raster_get_band(), rt_raster_get_geotransform_matrix(), rt_raster_get_height(), rt_raster_get_width(), rt_raster_new(), rt_raster_replace_band(), rt_raster_set_offsets(), rt_raster_set_scale(), rt_raster_set_skews(), rt_util_envelope_to_lwpoly(), rt_util_from_ogr_envelope(), rt_util_gdal_driver_registered(), rt_util_pixtype_to_gdal_datatype(), rtalloc(), rtdealloc(), rterror(), rtinfo(), rt_raster_t::scaleX, rt_raster_t::scaleY, _rti_rasterize_arg_t::src_sr, rt_envelope::UpperLeftX, rt_envelope::UpperLeftY, _rti_rasterize_arg_t::value, rt_raster_t::width, pixval::x, and pixval::y.

Referenced by RASTER_asRaster(), RASTER_clip(), RASTER_setPixelValuesGeomval(), and test_gdal_rasterize().

2521  {
2522  rt_raster rast = NULL;
2523  uint32_t i = 0;
2524  int err = 0;
2525 
2526  _rti_rasterize_arg arg = NULL;
2527 
2528  int _dim[2] = {0};
2529  double _scale[2] = {0};
2530  double _skew[2] = {0};
2531 
2532  OGRErr ogrerr;
2533  OGRGeometryH src_geom;
2534  OGREnvelope src_env;
2535  rt_envelope extent;
2536  OGRwkbGeometryType wkbtype = wkbUnknown;
2537 
2538  int ul_user = 0;
2539 
2540  CPLErr cplerr;
2541  double _gt[6] = {0};
2542  GDALDriverH _drv = NULL;
2543  int unload_drv = 0;
2544  GDALDatasetH _ds = NULL;
2545  GDALRasterBandH _band = NULL;
2546 
2547  uint16_t _width = 0;
2548  uint16_t _height = 0;
2549 
2550  RASTER_DEBUG(3, "starting");
2551 
2552  assert(NULL != wkb);
2553  assert(0 != wkb_len);
2554 
2555  /* internal variables */
2556  arg = _rti_rasterize_arg_init();
2557  if (arg == NULL) {
2558  rterror("rt_raster_gdal_rasterize: Could not initialize internal variables");
2559  return NULL;
2560  }
2561 
2562  /* no bands, raster is a mask */
2563  if (num_bands < 1) {
2564  arg->noband = 1;
2565  arg->numbands = 1;
2566 
2567  arg->pixtype = (rt_pixtype *) rtalloc(sizeof(rt_pixtype));
2568  arg->pixtype[0] = PT_8BUI;
2569 
2570  arg->init = (double *) rtalloc(sizeof(double));
2571  arg->init[0] = 0;
2572 
2573  arg->nodata = (double *) rtalloc(sizeof(double));
2574  arg->nodata[0] = 0;
2575 
2576  arg->hasnodata = (uint8_t *) rtalloc(sizeof(uint8_t));
2577  arg->hasnodata[0] = 1;
2578 
2579  arg->value = (double *) rtalloc(sizeof(double));
2580  arg->value[0] = 1;
2581  }
2582  else {
2583  arg->noband = 0;
2584  arg->numbands = num_bands;
2585 
2586  arg->pixtype = pixtype;
2587  arg->init = init;
2588  arg->nodata = nodata;
2589  arg->hasnodata = hasnodata;
2590  arg->value = value;
2591  }
2592 
2593  /* OGR spatial reference */
2594  if (NULL != srs && strlen(srs)) {
2595  arg->src_sr = OSRNewSpatialReference(NULL);
2596  if (OSRSetFromUserInput(arg->src_sr, srs) != OGRERR_NONE) {
2597  rterror("rt_raster_gdal_rasterize: Could not create OSR spatial reference using the provided srs: %s", srs);
2599  return NULL;
2600  }
2601  }
2602 
2603  /* convert WKB to OGR Geometry */
2604  ogrerr = OGR_G_CreateFromWkb((unsigned char *) wkb, arg->src_sr, &src_geom, wkb_len);
2605  if (ogrerr != OGRERR_NONE) {
2606  rterror("rt_raster_gdal_rasterize: Could not create OGR Geometry from WKB");
2607 
2609  /* OGRCleanupAll(); */
2610 
2611  return NULL;
2612  }
2613 
2614  /* OGR Geometry is empty */
2615  if (OGR_G_IsEmpty(src_geom)) {
2616  rtinfo("Geometry provided is empty. Returning empty raster");
2617 
2618  OGR_G_DestroyGeometry(src_geom);
2620  /* OGRCleanupAll(); */
2621 
2622  return rt_raster_new(0, 0);
2623  }
2624 
2625  /* get envelope */
2626  OGR_G_GetEnvelope(src_geom, &src_env);
2627  rt_util_from_ogr_envelope(src_env, &extent);
2628 
2629  RASTER_DEBUGF(3, "Suggested raster envelope: %f, %f, %f, %f",
2630  extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
2631 
2632  /* user-defined scale */
2633  if (
2634  (NULL != scale_x) &&
2635  (NULL != scale_y) &&
2636  (FLT_NEQ(*scale_x, 0.0)) &&
2637  (FLT_NEQ(*scale_y, 0.0))
2638  ) {
2639  /* for now, force scale to be in left-right, top-down orientation */
2640  _scale[0] = fabs(*scale_x);
2641  _scale[1] = fabs(*scale_y);
2642  }
2643  /* user-defined width/height */
2644  else if (
2645  (NULL != width) &&
2646  (NULL != height) &&
2647  (FLT_NEQ(*width, 0.0)) &&
2648  (FLT_NEQ(*height, 0.0))
2649  ) {
2650  _dim[0] = abs(*width);
2651  _dim[1] = abs(*height);
2652 
2653  if (FLT_NEQ(extent.MaxX, extent.MinX))
2654  _scale[0] = fabs((extent.MaxX - extent.MinX) / _dim[0]);
2655  else
2656  _scale[0] = 1.;
2657 
2658  if (FLT_NEQ(extent.MaxY, extent.MinY))
2659  _scale[1] = fabs((extent.MaxY - extent.MinY) / _dim[1]);
2660  else
2661  _scale[1] = 1.;
2662  }
2663  else {
2664  rterror("rt_raster_gdal_rasterize: Values must be provided for width and height or X and Y of scale");
2665 
2666  OGR_G_DestroyGeometry(src_geom);
2668  /* OGRCleanupAll(); */
2669 
2670  return NULL;
2671  }
2672  RASTER_DEBUGF(3, "scale (x, y) = %f, %f", _scale[0], -1 * _scale[1]);
2673  RASTER_DEBUGF(3, "dim (x, y) = %d, %d", _dim[0], _dim[1]);
2674 
2675  /* user-defined skew */
2676  if (NULL != skew_x) {
2677  _skew[0] = *skew_x;
2678 
2679  /*
2680  negative scale-x affects skew
2681  for now, force skew to be in left-right, top-down orientation
2682  */
2683  if (
2684  NULL != scale_x &&
2685  *scale_x < 0.
2686  ) {
2687  _skew[0] *= -1;
2688  }
2689  }
2690  if (NULL != skew_y) {
2691  _skew[1] = *skew_y;
2692 
2693  /*
2694  positive scale-y affects skew
2695  for now, force skew to be in left-right, top-down orientation
2696  */
2697  if (
2698  NULL != scale_y &&
2699  *scale_y > 0.
2700  ) {
2701  _skew[1] *= -1;
2702  }
2703  }
2704 
2705  /*
2706  if geometry is a point, a linestring or set of either and bounds not set,
2707  increase extent by a pixel to avoid missing points on border
2708 
2709  a whole pixel is used instead of half-pixel due to backward
2710  compatibility with GDAL 1.6, 1.7 and 1.8. 1.9+ works fine with half-pixel.
2711  */
2712  wkbtype = wkbFlatten(OGR_G_GetGeometryType(src_geom));
2713  if ((
2714  (wkbtype == wkbPoint) ||
2715  (wkbtype == wkbMultiPoint) ||
2716  (wkbtype == wkbLineString) ||
2717  (wkbtype == wkbMultiLineString)
2718  ) &&
2719  _dim[0] == 0 &&
2720  _dim[1] == 0
2721  ) {
2722  int result;
2723  LWPOLY *epoly = NULL;
2724  LWGEOM *lwgeom = NULL;
2725  GEOSGeometry *egeom = NULL;
2726  GEOSGeometry *geom = NULL;
2727 
2728  RASTER_DEBUG(3, "Testing geometry is properly contained by extent");
2729 
2730  /*
2731  see if geometry is properly contained by extent
2732  all parts of geometry lies within extent
2733  */
2734 
2735  /* initialize GEOS */
2736  initGEOS(rtinfo, lwgeom_geos_error);
2737 
2738  /* convert envelope to geometry */
2739  RASTER_DEBUG(4, "Converting envelope to geometry");
2740  epoly = rt_util_envelope_to_lwpoly(extent);
2741  if (epoly == NULL) {
2742  rterror("rt_raster_gdal_rasterize: Could not create envelope's geometry to test if geometry is properly contained by extent");
2743 
2744  OGR_G_DestroyGeometry(src_geom);
2746  /* OGRCleanupAll(); */
2747 
2748  return NULL;
2749  }
2750 
2751  egeom = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(epoly), 0);
2752  lwpoly_free(epoly);
2753 
2754  /* convert WKB to geometry */
2755  RASTER_DEBUG(4, "Converting WKB to geometry");
2756  lwgeom = lwgeom_from_wkb(wkb, wkb_len, LW_PARSER_CHECK_NONE);
2757  geom = (GEOSGeometry *) LWGEOM2GEOS(lwgeom, 0);
2758  lwgeom_free(lwgeom);
2759 
2760  result = GEOSRelatePattern(egeom, geom, "T**FF*FF*");
2761  GEOSGeom_destroy(geom);
2762  GEOSGeom_destroy(egeom);
2763 
2764  if (result == 2) {
2765  rterror("rt_raster_gdal_rasterize: Could not test if geometry is properly contained by extent for geometry within extent");
2766 
2767  OGR_G_DestroyGeometry(src_geom);
2769  /* OGRCleanupAll(); */
2770 
2771  return NULL;
2772  }
2773 
2774  /* geometry NOT properly contained by extent */
2775  if (!result) {
2776 
2777 #if POSTGIS_GDAL_VERSION > 18
2778 
2779  /* check alignment flag: grid_xw */
2780  if (
2781  (NULL == ul_xw && NULL == ul_yw) &&
2782  (NULL != grid_xw && NULL != grid_yw) &&
2783  FLT_NEQ(*grid_xw, extent.MinX)
2784  ) {
2785  /* do nothing */
2786  RASTER_DEBUG(3, "Skipping extent adjustment on X-axis due to upcoming alignment");
2787  }
2788  else {
2789  RASTER_DEBUG(3, "Adjusting extent for GDAL > 1.8 by half the scale on X-axis");
2790  extent.MinX -= (_scale[0] / 2.);
2791  extent.MaxX += (_scale[0] / 2.);
2792  }
2793 
2794  /* check alignment flag: grid_yw */
2795  if (
2796  (NULL == ul_xw && NULL == ul_yw) &&
2797  (NULL != grid_xw && NULL != grid_yw) &&
2798  FLT_NEQ(*grid_yw, extent.MaxY)
2799  ) {
2800  /* do nothing */
2801  RASTER_DEBUG(3, "Skipping extent adjustment on Y-axis due to upcoming alignment");
2802  }
2803  else {
2804  RASTER_DEBUG(3, "Adjusting extent for GDAL > 1.8 by half the scale on Y-axis");
2805  extent.MinY -= (_scale[1] / 2.);
2806  extent.MaxY += (_scale[1] / 2.);
2807  }
2808 
2809 #else
2810 
2811  /* check alignment flag: grid_xw */
2812  if (
2813  (NULL == ul_xw && NULL == ul_yw) &&
2814  (NULL != grid_xw && NULL != grid_yw) &&
2815  FLT_NEQ(*grid_xw, extent.MinX)
2816  ) {
2817  /* do nothing */
2818  RASTER_DEBUG(3, "Skipping extent adjustment on X-axis due to upcoming alignment");
2819  }
2820  else {
2821  RASTER_DEBUG(3, "Adjusting extent for GDAL <= 1.8 by the scale on X-axis");
2822  extent.MinX -= _scale[0];
2823  extent.MaxX += _scale[0];
2824  }
2825 
2826 
2827  /* check alignment flag: grid_yw */
2828  if (
2829  (NULL == ul_xw && NULL == ul_yw) &&
2830  (NULL != grid_xw && NULL != grid_yw) &&
2831  FLT_NEQ(*grid_yw, extent.MaxY)
2832  ) {
2833  /* do nothing */
2834  RASTER_DEBUG(3, "Skipping extent adjustment on Y-axis due to upcoming alignment");
2835  }
2836  else {
2837  RASTER_DEBUG(3, "Adjusting extent for GDAL <= 1.8 by the scale on Y-axis");
2838  extent.MinY -= _scale[1];
2839  extent.MaxY += _scale[1];
2840  }
2841 
2842 #endif
2843 
2844  }
2845 
2846  RASTER_DEBUGF(3, "Adjusted extent: %f, %f, %f, %f",
2847  extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
2848 
2849  extent.UpperLeftX = extent.MinX;
2850  extent.UpperLeftY = extent.MaxY;
2851  }
2852 
2853  /* reprocess extent if skewed */
2854  if (
2855  FLT_NEQ(_skew[0], 0) ||
2856  FLT_NEQ(_skew[1], 0)
2857  ) {
2858  rt_raster skewedrast;
2859 
2860  RASTER_DEBUG(3, "Computing skewed extent's envelope");
2861 
2862  skewedrast = rt_raster_compute_skewed_raster(
2863  extent,
2864  _skew,
2865  _scale,
2866  0.01
2867  );
2868  if (skewedrast == NULL) {
2869  rterror("rt_raster_gdal_rasterize: Could not compute skewed raster");
2870 
2871  OGR_G_DestroyGeometry(src_geom);
2873  /* OGRCleanupAll(); */
2874 
2875  return NULL;
2876  }
2877 
2878  _dim[0] = skewedrast->width;
2879  _dim[1] = skewedrast->height;
2880 
2881  extent.UpperLeftX = skewedrast->ipX;
2882  extent.UpperLeftY = skewedrast->ipY;
2883 
2884  rt_raster_destroy(skewedrast);
2885  }
2886 
2887  /* raster dimensions */
2888  if (!_dim[0])
2889  _dim[0] = (int) fmax((fabs(extent.MaxX - extent.MinX) + (_scale[0] / 2.)) / _scale[0], 1);
2890  if (!_dim[1])
2891  _dim[1] = (int) fmax((fabs(extent.MaxY - extent.MinY) + (_scale[1] / 2.)) / _scale[1], 1);
2892 
2893  /* temporary raster */
2894  rast = rt_raster_new(_dim[0], _dim[1]);
2895  if (rast == NULL) {
2896  rterror("rt_raster_gdal_rasterize: Out of memory allocating temporary raster");
2897 
2898  OGR_G_DestroyGeometry(src_geom);
2900  /* OGRCleanupAll(); */
2901 
2902  return NULL;
2903  }
2904 
2905  /* set raster's spatial attributes */
2906  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
2907  rt_raster_set_scale(rast, _scale[0], -1 * _scale[1]);
2908  rt_raster_set_skews(rast, _skew[0], _skew[1]);
2909 
2911  RASTER_DEBUGF(3, "Temp raster's geotransform: %f, %f, %f, %f, %f, %f",
2912  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
2913  RASTER_DEBUGF(3, "Temp raster's dimensions (width x height): %d x %d",
2914  _dim[0], _dim[1]);
2915 
2916  /* user-specified upper-left corner */
2917  if (
2918  NULL != ul_xw &&
2919  NULL != ul_yw
2920  ) {
2921  ul_user = 1;
2922 
2923  RASTER_DEBUGF(4, "Using user-specified upper-left corner: %f, %f", *ul_xw, *ul_yw);
2924 
2925  /* set upper-left corner */
2926  rt_raster_set_offsets(rast, *ul_xw, *ul_yw);
2927  extent.UpperLeftX = *ul_xw;
2928  extent.UpperLeftY = *ul_yw;
2929  }
2930  else if (
2931  ((NULL != ul_xw) && (NULL == ul_yw)) ||
2932  ((NULL == ul_xw) && (NULL != ul_yw))
2933  ) {
2934  rterror("rt_raster_gdal_rasterize: Both X and Y upper-left corner values must be provided");
2935 
2936  rt_raster_destroy(rast);
2937  OGR_G_DestroyGeometry(src_geom);
2939  /* OGRCleanupAll(); */
2940 
2941  return NULL;
2942  }
2943 
2944  /* alignment only considered if upper-left corner not provided */
2945  if (
2946  !ul_user && (
2947  (NULL != grid_xw) || (NULL != grid_yw)
2948  )
2949  ) {
2950 
2951  if (
2952  ((NULL != grid_xw) && (NULL == grid_yw)) ||
2953  ((NULL == grid_xw) && (NULL != grid_yw))
2954  ) {
2955  rterror("rt_raster_gdal_rasterize: Both X and Y alignment values must be provided");
2956 
2957  rt_raster_destroy(rast);
2958  OGR_G_DestroyGeometry(src_geom);
2960  /* OGRCleanupAll(); */
2961 
2962  return NULL;
2963  }
2964 
2965  RASTER_DEBUGF(4, "Aligning extent to user-specified grid: %f, %f", *grid_xw, *grid_yw);
2966 
2967  do {
2968  double _r[2] = {0};
2969  double _w[2] = {0};
2970 
2971  /* raster is already aligned */
2972  if (FLT_EQ(*grid_xw, extent.UpperLeftX) && FLT_EQ(*grid_yw, extent.UpperLeftY)) {
2973  RASTER_DEBUG(3, "Skipping raster alignment as it is already aligned to grid");
2974  break;
2975  }
2976 
2977  extent.UpperLeftX = rast->ipX;
2978  extent.UpperLeftY = rast->ipY;
2979  rt_raster_set_offsets(rast, *grid_xw, *grid_yw);
2980 
2981  /* process upper-left corner */
2983  rast,
2984  extent.UpperLeftX, extent.UpperLeftY,
2985  &(_r[0]), &(_r[1]),
2986  NULL
2987  ) != ES_NONE) {
2988  rterror("rt_raster_gdal_rasterize: Could not compute raster pixel for spatial coordinates");
2989 
2990  rt_raster_destroy(rast);
2991  OGR_G_DestroyGeometry(src_geom);
2993  /* OGRCleanupAll(); */
2994 
2995  return NULL;
2996  }
2997 
2999  rast,
3000  _r[0], _r[1],
3001  &(_w[0]), &(_w[1]),
3002  NULL
3003  ) != ES_NONE) {
3004  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3005 
3006  rt_raster_destroy(rast);
3007  OGR_G_DestroyGeometry(src_geom);
3009  /* OGRCleanupAll(); */
3010 
3011  return NULL;
3012  }
3013 
3014  /* shift occurred */
3015  if (FLT_NEQ(_w[0], extent.UpperLeftX)) {
3016  if (NULL == width)
3017  rast->width++;
3018  else if (NULL == scale_x) {
3019  double _c[2] = {0};
3020 
3021  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
3022 
3023  /* get upper-right corner */
3025  rast,
3026  rast->width, 0,
3027  &(_c[0]), &(_c[1]),
3028  NULL
3029  ) != ES_NONE) {
3030  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3031 
3032  rt_raster_destroy(rast);
3033  OGR_G_DestroyGeometry(src_geom);
3035  /* OGRCleanupAll(); */
3036 
3037  return NULL;
3038  }
3039 
3040  rast->scaleX = fabs((_c[0] - _w[0]) / ((double) rast->width));
3041  }
3042  }
3043  if (FLT_NEQ(_w[1], extent.UpperLeftY)) {
3044  if (NULL == height)
3045  rast->height++;
3046  else if (NULL == scale_y) {
3047  double _c[2] = {0};
3048 
3049  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
3050 
3051  /* get upper-right corner */
3053  rast,
3054  0, rast->height,
3055  &(_c[0]), &(_c[1]),
3056  NULL
3057  ) != ES_NONE) {
3058  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3059 
3060  rt_raster_destroy(rast);
3061  OGR_G_DestroyGeometry(src_geom);
3063  /* OGRCleanupAll(); */
3064 
3065  return NULL;
3066  }
3067 
3068  rast->scaleY = -1 * fabs((_c[1] - _w[1]) / ((double) rast->height));
3069  }
3070  }
3071 
3072  rt_raster_set_offsets(rast, _w[0], _w[1]);
3073  }
3074  while (0);
3075  }
3076 
3077  /*
3078  after this point, rt_envelope extent is no longer used
3079  */
3080 
3081  /* get key attributes from rast */
3082  _dim[0] = rast->width;
3083  _dim[1] = rast->height;
3085 
3086  /* scale-x is negative or scale-y is positive */
3087  if ((
3088  (NULL != scale_x) && (*scale_x < 0.)
3089  ) || (
3090  (NULL != scale_y) && (*scale_y > 0)
3091  )) {
3092  double _w[2] = {0};
3093 
3094  /* negative scale-x */
3095  if (
3096  (NULL != scale_x) &&
3097  (*scale_x < 0.)
3098  ) {
3099  RASTER_DEBUG(3, "Processing negative scale-x");
3100 
3102  rast,
3103  _dim[0], 0,
3104  &(_w[0]), &(_w[1]),
3105  NULL
3106  ) != ES_NONE) {
3107  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3108 
3109  rt_raster_destroy(rast);
3110  OGR_G_DestroyGeometry(src_geom);
3112  /* OGRCleanupAll(); */
3113 
3114  return NULL;
3115  }
3116 
3117  _gt[0] = _w[0];
3118  _gt[1] = *scale_x;
3119 
3120  /* check for skew */
3121  if (NULL != skew_x && FLT_NEQ(*skew_x, 0))
3122  _gt[2] = *skew_x;
3123  }
3124  /* positive scale-y */
3125  if (
3126  (NULL != scale_y) &&
3127  (*scale_y > 0)
3128  ) {
3129  RASTER_DEBUG(3, "Processing positive scale-y");
3130 
3132  rast,
3133  0, _dim[1],
3134  &(_w[0]), &(_w[1]),
3135  NULL
3136  ) != ES_NONE) {
3137  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
3138 
3139  rt_raster_destroy(rast);
3140  OGR_G_DestroyGeometry(src_geom);
3142  /* OGRCleanupAll(); */
3143 
3144  return NULL;
3145  }
3146 
3147  _gt[3] = _w[1];
3148  _gt[5] = *scale_y;
3149 
3150  /* check for skew */
3151  if (NULL != skew_y && FLT_NEQ(*skew_y, 0))
3152  _gt[4] = *skew_y;
3153  }
3154  }
3155 
3156  rt_raster_destroy(rast);
3157  rast = NULL;
3158 
3159  RASTER_DEBUGF(3, "Applied geotransform: %f, %f, %f, %f, %f, %f",
3160  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
3161  RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d",
3162  _dim[0], _dim[1]);
3163 
3164  /* load GDAL mem */
3165  if (!rt_util_gdal_driver_registered("MEM")) {
3166  RASTER_DEBUG(4, "Registering MEM driver");
3167  GDALRegister_MEM();
3168  unload_drv = 1;
3169  }
3170  _drv = GDALGetDriverByName("MEM");
3171  if (NULL == _drv) {
3172  rterror("rt_raster_gdal_rasterize: Could not load the MEM GDAL driver");
3173 
3174  OGR_G_DestroyGeometry(src_geom);
3176  /* OGRCleanupAll(); */
3177 
3178  return NULL;
3179  }
3180 
3181  /* unload driver from GDAL driver manager */
3182  if (unload_drv) {
3183  RASTER_DEBUG(4, "Deregistering MEM driver");
3184  GDALDeregisterDriver(_drv);
3185  }
3186 
3187  _ds = GDALCreate(_drv, "", _dim[0], _dim[1], 0, GDT_Byte, NULL);
3188  if (NULL == _ds) {
3189  rterror("rt_raster_gdal_rasterize: Could not create a GDALDataset to rasterize the geometry into");
3190 
3191  OGR_G_DestroyGeometry(src_geom);
3193  /* OGRCleanupAll(); */
3194  if (unload_drv) GDALDestroyDriver(_drv);
3195 
3196  return NULL;
3197  }
3198 
3199  /* set geotransform */
3200  cplerr = GDALSetGeoTransform(_ds, _gt);
3201  if (cplerr != CE_None) {
3202  rterror("rt_raster_gdal_rasterize: Could not set geotransform on GDALDataset");
3203 
3204  OGR_G_DestroyGeometry(src_geom);
3206  /* OGRCleanupAll(); */
3207 
3208  GDALClose(_ds);
3209  if (unload_drv) GDALDestroyDriver(_drv);
3210 
3211  return NULL;
3212  }
3213 
3214  /* set SRS */
3215  if (NULL != arg->src_sr) {
3216  char *_srs = NULL;
3217  OSRExportToWkt(arg->src_sr, &_srs);
3218 
3219  cplerr = GDALSetProjection(_ds, _srs);
3220  CPLFree(_srs);
3221  if (cplerr != CE_None) {
3222  rterror("rt_raster_gdal_rasterize: Could not set projection on GDALDataset");
3223 
3224  OGR_G_DestroyGeometry(src_geom);
3226  /* OGRCleanupAll(); */
3227 
3228  GDALClose(_ds);
3229  if (unload_drv) GDALDestroyDriver(_drv);
3230 
3231  return NULL;
3232  }
3233  }
3234 
3235  /* set bands */
3236  for (i = 0; i < arg->numbands; i++) {
3237  err = 0;
3238 
3239  do {
3240  /* add band */
3241  cplerr = GDALAddBand(_ds, rt_util_pixtype_to_gdal_datatype(arg->pixtype[i]), NULL);
3242  if (cplerr != CE_None) {
3243  rterror("rt_raster_gdal_rasterize: Could not add band to GDALDataset");
3244  err = 1;
3245  break;
3246  }
3247 
3248  _band = GDALGetRasterBand(_ds, i + 1);
3249  if (NULL == _band) {
3250  rterror("rt_raster_gdal_rasterize: Could not get band %d from GDALDataset", i + 1);
3251  err = 1;
3252  break;
3253  }
3254 
3255  /* nodata value */
3256  if (arg->hasnodata[i]) {
3257  RASTER_DEBUGF(4, "Setting NODATA value of band %d to %f", i, arg->nodata[i]);
3258  cplerr = GDALSetRasterNoDataValue(_band, arg->nodata[i]);
3259  if (cplerr != CE_None) {
3260  rterror("rt_raster_gdal_rasterize: Could not set nodata value");
3261  err = 1;
3262  break;
3263  }
3264  RASTER_DEBUGF(4, "NODATA value set to %f", GDALGetRasterNoDataValue(_band, NULL));
3265  }
3266 
3267  /* initial value */
3268  cplerr = GDALFillRaster(_band, arg->init[i], 0);
3269  if (cplerr != CE_None) {
3270  rterror("rt_raster_gdal_rasterize: Could not set initial value");
3271  err = 1;
3272  break;
3273  }
3274  }
3275  while (0);
3276 
3277  if (err) {
3278 
3279  OGR_G_DestroyGeometry(src_geom);
3281 
3282  /* OGRCleanupAll(); */
3283 
3284  GDALClose(_ds);
3285  if (unload_drv) GDALDestroyDriver(_drv);
3286 
3287  return NULL;
3288  }
3289  }
3290 
3291  arg->bandlist = (int *) rtalloc(sizeof(int) * arg->numbands);
3292  for (i = 0; i < arg->numbands; i++) arg->bandlist[i] = i + 1;
3293 
3294  /* burn geometry */
3295  cplerr = GDALRasterizeGeometries(
3296  _ds,
3297  arg->numbands, arg->bandlist,
3298  1, &src_geom,
3299  NULL, NULL,
3300  arg->value,
3301  options,
3302  NULL, NULL
3303  );
3304  if (cplerr != CE_None) {
3305  rterror("rt_raster_gdal_rasterize: Could not rasterize geometry");
3306 
3307  OGR_G_DestroyGeometry(src_geom);
3309  /* OGRCleanupAll(); */
3310 
3311  GDALClose(_ds);
3312  if (unload_drv) GDALDestroyDriver(_drv);
3313 
3314  return NULL;
3315  }
3316 
3317  /* convert gdal dataset to raster */
3318  GDALFlushCache(_ds);
3319  RASTER_DEBUG(3, "Converting GDAL dataset to raster");
3320  rast = rt_raster_from_gdal_dataset(_ds);
3321 
3322  OGR_G_DestroyGeometry(src_geom);
3323  /* OGRCleanupAll(); */
3324 
3325  GDALClose(_ds);
3326  if (unload_drv) GDALDestroyDriver(_drv);
3327 
3328  if (NULL == rast) {
3329  rterror("rt_raster_gdal_rasterize: Could not rasterize geometry");
3330  return NULL;
3331  }
3332 
3333  /* width, height */
3334  _width = rt_raster_get_width(rast);
3335  _height = rt_raster_get_height(rast);
3336 
3337  /* check each band for pixtype */
3338  for (i = 0; i < arg->numbands; i++) {
3339  uint8_t *data = NULL;
3340  rt_band band = NULL;
3341  rt_band oldband = NULL;
3342 
3343  double val = 0;
3344  int nodata = 0;
3345  int hasnodata = 0;
3346  double nodataval = 0;
3347  int x = 0;
3348  int y = 0;
3349 
3350  oldband = rt_raster_get_band(rast, i);
3351  if (oldband == NULL) {
3352  rterror("rt_raster_gdal_rasterize: Could not get band %d of output raster", i);
3354  rt_raster_destroy(rast);
3355  return NULL;
3356  }
3357 
3358  /* band is of user-specified type */
3359  if (rt_band_get_pixtype(oldband) == arg->pixtype[i])
3360  continue;
3361 
3362  /* hasnodata, nodataval */
3363  hasnodata = rt_band_get_hasnodata_flag(oldband);
3364  if (hasnodata)
3365  rt_band_get_nodata(oldband, &nodataval);
3366 
3367  /* allocate data */
3368  data = rtalloc(rt_pixtype_size(arg->pixtype[i]) * _width * _height);
3369  if (data == NULL) {
3370  rterror("rt_raster_gdal_rasterize: Could not allocate memory for band data");
3372  rt_raster_destroy(rast);
3373  return NULL;
3374  }
3375  memset(data, 0, rt_pixtype_size(arg->pixtype[i]) * _width * _height);
3376 
3377  /* create new band of correct type */
3378  band = rt_band_new_inline(
3379  _width, _height,
3380  arg->pixtype[i],
3381  hasnodata, nodataval,
3382  data
3383  );
3384  if (band == NULL) {
3385  rterror("rt_raster_gdal_rasterize: Could not create band");
3386  rtdealloc(data);
3388  rt_raster_destroy(rast);
3389  return NULL;
3390  }
3391 
3392  /* give ownership of data to band */
3393  rt_band_set_ownsdata_flag(band, 1);
3394 
3395  /* copy pixel by pixel */
3396  for (x = 0; x < _width; x++) {
3397  for (y = 0; y < _height; y++) {
3398  err = rt_band_get_pixel(oldband, x, y, &val, &nodata);
3399  if (err != ES_NONE) {
3400  rterror("rt_raster_gdal_rasterize: Could not get pixel value");
3402  rt_raster_destroy(rast);
3403  rt_band_destroy(band);
3404  return NULL;
3405  }
3406 
3407  if (nodata)
3408  val = nodataval;
3409 
3410  err = rt_band_set_pixel(band, x, y, val, NULL);
3411  if (err != ES_NONE) {
3412  rterror("rt_raster_gdal_rasterize: Could not set pixel value");
3414  rt_raster_destroy(rast);
3415  rt_band_destroy(band);
3416  return NULL;
3417  }
3418  }
3419  }
3420 
3421  /* replace band */
3422  oldband = rt_raster_replace_band(rast, band, i);
3423  if (oldband == NULL) {
3424  rterror("rt_raster_gdal_rasterize: Could not replace band %d of output raster", i);
3426  rt_raster_destroy(rast);
3427  rt_band_destroy(band);
3428  return NULL;
3429  }
3430 
3431  /* free oldband */
3432  rt_band_destroy(oldband);
3433  }
3434 
3436 
3437  RASTER_DEBUG(3, "done");
3438 
3439  return rast;
3440 }
double MinY
Definition: librtcore.h:167
double UpperLeftX
Definition: librtcore.h:170
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_raster.c:137
rt_band rt_raster_replace_band(rt_raster raster, rt_band band, int index)
Replace band at provided index with new band.
Definition: rt_raster.c:1498
double MaxY
Definition: librtcore.h:168
double MaxX
Definition: librtcore.h:166
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1144
band
Definition: ovdump.py:57
rt_errorstate rt_raster_cell_to_geopoint(rt_raster raster, double xr, double yr, double *xw, double *yw, double *gt)
Convert an xr, yr raster point to an xw, yw point on map.
Definition: rt_raster.c:755
#define FLT_EQ(x, y)
Definition: librtcore.h:2234
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
data
Definition: ovdump.py:103
LWPOLY * rt_util_envelope_to_lwpoly(rt_envelope ext)
Definition: rt_util.c:437
OGRSpatialReferenceH src_sr
Definition: rt_raster.c:2422
double ipY
Definition: librtcore.h:2296
rt_errorstate rt_raster_geopoint_to_cell(rt_raster raster, double xw, double yw, double *xr, double *yr, double *igt)
Convert an xw,yw map point to a xr,yr raster point.
Definition: rt_raster.c:806
rt_raster rt_raster_compute_skewed_raster(rt_envelope extent, double *skew, double *scale, double tolerance)
Definition: rt_raster.c:960
rt_pixtype
Definition: librtcore.h:185
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:340
static _rti_rasterize_arg _rti_rasterize_arg_init()
Definition: rt_raster.c:2433
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:320
uint16_t height
Definition: librtcore.h:2302
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1730
unsigned int uint32_t
Definition: uthash.h:78
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:2004
rt_raster rt_raster_from_gdal_dataset(GDALDatasetH ds)
Return a raster from a GDAL dataset.
Definition: rt_raster.c:2182
static void _rti_rasterize_arg_destroy(_rti_rasterize_arg arg)
Definition: rt_raster.c:2459
void rt_band_set_ownsdata_flag(rt_band band, int flag)
Definition: rt_band.c:667
GDALDataType rt_util_pixtype_to_gdal_datatype(rt_pixtype pt)
Convert rt_pixtype to GDALDataType.
Definition: rt_util.c:118
int rt_util_gdal_driver_registered(const char *drv)
Definition: rt_util.c:355
void lwgeom_geos_error(const char *fmt,...)
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:175
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1221
uint16_t width
Definition: librtcore.h:2301
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_raster.c:706
#define FLT_NEQ(x, y)
Definition: librtcore.h:2233
void rtinfo(const char *fmt,...)
Definition: rt_context.c:211
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:674
rt_pixtype * pixtype
Definition: rt_raster.c:2424
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_raster.c:199
rt_band rt_band_new_inline(uint16_t width, uint16_t height, rt_pixtype pixtype, uint32_t hasnodata, double nodataval, uint8_t *data)
Create an in-db rt_band with no data.
Definition: rt_band.c:63
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, uint8_t autofix)
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition: rt_pixel.c:39
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
void rtdealloc(void *mem)
Definition: rt_context.c:186
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
Definition: rt_raster.c:168
double MinX
Definition: librtcore.h:165
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:631
void rt_util_from_ogr_envelope(OGREnvelope env, rt_envelope *ext)
Definition: rt_util.c:408
double ipX
Definition: librtcore.h:2295
rt_errorstate rt_band_set_pixel(rt_band band, int x, int y, double val, int *converted)
Set single pixel&#39;s value.
Definition: rt_band.c:974
int value
Definition: genraster.py:61
double scaleX
Definition: librtcore.h:2293
LWGEOM * lwgeom_from_wkb(const uint8_t *wkb, const size_t wkb_size, const char check)
WKB inputs must have a declared size, to prevent malformed WKB from reading off the end of the memory...
Definition: lwin_wkb.c:770
unsigned char uint8_t
Definition: uthash.h:79
double UpperLeftY
Definition: librtcore.h:171
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
double scaleY
Definition: librtcore.h:2294
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
Here is the call graph for this function:
Here is the caller graph for this function: