PostGIS  3.0.6dev-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 2504 of file rt_raster.c.

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

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, 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_from_ogr_envelope(), rt_util_gdal_driver_registered(), rt_util_pixtype_to_gdal_datatype(), rtalloc(), rtdealloc(), rterror(), rtinfo(), _rti_rasterize_arg_t::src_sr, rt_envelope::UpperLeftX, rt_envelope::UpperLeftY, _rti_rasterize_arg_t::value, genraster::value, rt_raster_t::width, pixval::x, and pixval::y.

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

Here is the call graph for this function:
Here is the caller graph for this function: