PostGIS  3.1.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 2491 of file rt_raster.c.

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