PostGIS  2.5.0beta1dev-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 2513 of file rt_raster.c.

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

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

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