PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ rt_raster_gdal_rasterize()

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

Return a raster of the provided geometry.

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

Definition at line 2509 of file rt_raster.c.

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

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: