PostGIS  2.4.9dev-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 2492 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().

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