PostGIS  2.1.10dev-r@@SVN_REVISION@@
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 10620 of file rt_api.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(), lwnotice(), 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, result, 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, rt_envelope::UpperLeftX, rt_envelope::UpperLeftY, genraster::value, _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().

10632  {
10633  rt_raster rast = NULL;
10634  int i = 0;
10635  int err = 0;
10636 
10637  _rti_rasterize_arg arg = NULL;
10638 
10639  int _dim[2] = {0};
10640  double _scale[2] = {0};
10641  double _skew[2] = {0};
10642 
10643  OGRErr ogrerr;
10644  OGRSpatialReferenceH src_sr = NULL;
10645  OGRGeometryH src_geom;
10646  OGREnvelope src_env;
10647  rt_envelope extent;
10648  OGRwkbGeometryType wkbtype = wkbUnknown;
10649 
10650  int ul_user = 0;
10651 
10652  CPLErr cplerr;
10653  double _gt[6] = {0};
10654  GDALDriverH _drv = NULL;
10655  int unload_drv = 0;
10656  GDALDatasetH _ds = NULL;
10657  GDALRasterBandH _band = NULL;
10658 
10659  uint16_t _width = 0;
10660  uint16_t _height = 0;
10661 
10662  RASTER_DEBUG(3, "starting");
10663 
10664  assert(NULL != wkb);
10665  assert(0 != wkb_len);
10666 
10667  /* internal variables */
10668  arg = _rti_rasterize_arg_init();
10669  if (arg == NULL) {
10670  rterror("rt_raster_gdal_rasterize: Could not initialize internal variables");
10671  return NULL;
10672  }
10673 
10674  /* no bands, raster is a mask */
10675  if (num_bands < 1) {
10676  arg->noband = 1;
10677  arg->numbands = 1;
10678 
10679  arg->pixtype = (rt_pixtype *) rtalloc(sizeof(rt_pixtype));
10680  arg->pixtype[0] = PT_8BUI;
10681 
10682  arg->init = (double *) rtalloc(sizeof(double));
10683  arg->init[0] = 0;
10684 
10685  arg->nodata = (double *) rtalloc(sizeof(double));
10686  arg->nodata[0] = 0;
10687 
10688  arg->hasnodata = (uint8_t *) rtalloc(sizeof(uint8_t));
10689  arg->hasnodata[0] = 1;
10690 
10691  arg->value = (double *) rtalloc(sizeof(double));
10692  arg->value[0] = 1;
10693  }
10694  else {
10695  arg->noband = 0;
10696  arg->numbands = num_bands;
10697 
10698  arg->pixtype = pixtype;
10699  arg->init = init;
10700  arg->nodata = nodata;
10701  arg->hasnodata = hasnodata;
10702  arg->value = value;
10703  }
10704 
10705  /* OGR spatial reference */
10706  if (NULL != srs && strlen(srs)) {
10707  src_sr = OSRNewSpatialReference(NULL);
10708  if (OSRSetFromUserInput(src_sr, srs) != OGRERR_NONE) {
10709  rterror("rt_raster_gdal_rasterize: Could not create OSR spatial reference using the provided srs: %s", srs);
10711  return NULL;
10712  }
10713  }
10714 
10715  /* convert WKB to OGR Geometry */
10716  ogrerr = OGR_G_CreateFromWkb((unsigned char *) wkb, src_sr, &src_geom, wkb_len);
10717  if (ogrerr != OGRERR_NONE) {
10718  rterror("rt_raster_gdal_rasterize: Could not create OGR Geometry from WKB");
10719 
10721 
10722  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
10723  /* OGRCleanupAll(); */
10724 
10725  return NULL;
10726  }
10727 
10728  /* OGR Geometry is empty */
10729  if (OGR_G_IsEmpty(src_geom)) {
10730  rtinfo("Geometry provided is empty. Returning empty raster");
10731 
10733 
10734  OGR_G_DestroyGeometry(src_geom);
10735  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
10736  /* OGRCleanupAll(); */
10737 
10738  return rt_raster_new(0, 0);
10739  }
10740 
10741  /* get envelope */
10742  OGR_G_GetEnvelope(src_geom, &src_env);
10743  rt_util_from_ogr_envelope(src_env, &extent);
10744 
10745  RASTER_DEBUGF(3, "Suggested raster envelope: %f, %f, %f, %f",
10746  extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
10747 
10748  /* user-defined scale */
10749  if (
10750  (NULL != scale_x) &&
10751  (NULL != scale_y) &&
10752  (FLT_NEQ(*scale_x, 0.0)) &&
10753  (FLT_NEQ(*scale_y, 0.0))
10754  ) {
10755  /* for now, force scale to be in left-right, top-down orientation */
10756  _scale[0] = fabs(*scale_x);
10757  _scale[1] = fabs(*scale_y);
10758  }
10759  /* user-defined width/height */
10760  else if (
10761  (NULL != width) &&
10762  (NULL != height) &&
10763  (FLT_NEQ(*width, 0.0)) &&
10764  (FLT_NEQ(*height, 0.0))
10765  ) {
10766  _dim[0] = abs(*width);
10767  _dim[1] = abs(*height);
10768 
10769  if (FLT_NEQ(extent.MaxX, extent.MinX))
10770  _scale[0] = fabs((extent.MaxX - extent.MinX) / _dim[0]);
10771  else
10772  _scale[0] = 1.;
10773 
10774  if (FLT_NEQ(extent.MaxY, extent.MinY))
10775  _scale[1] = fabs((extent.MaxY - extent.MinY) / _dim[1]);
10776  else
10777  _scale[1] = 1.;
10778  }
10779  else {
10780  rterror("rt_raster_gdal_rasterize: Values must be provided for width and height or X and Y of scale");
10781 
10783 
10784  OGR_G_DestroyGeometry(src_geom);
10785  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
10786  /* OGRCleanupAll(); */
10787 
10788  return NULL;
10789  }
10790  RASTER_DEBUGF(3, "scale (x, y) = %f, %f", _scale[0], -1 * _scale[1]);
10791  RASTER_DEBUGF(3, "dim (x, y) = %d, %d", _dim[0], _dim[1]);
10792 
10793  /* user-defined skew */
10794  if (NULL != skew_x) {
10795  _skew[0] = *skew_x;
10796 
10797  /*
10798  negative scale-x affects skew
10799  for now, force skew to be in left-right, top-down orientation
10800  */
10801  if (
10802  NULL != scale_x &&
10803  *scale_x < 0.
10804  ) {
10805  _skew[0] *= -1;
10806  }
10807  }
10808  if (NULL != skew_y) {
10809  _skew[1] = *skew_y;
10810 
10811  /*
10812  positive scale-y affects skew
10813  for now, force skew to be in left-right, top-down orientation
10814  */
10815  if (
10816  NULL != scale_y &&
10817  *scale_y > 0.
10818  ) {
10819  _skew[1] *= -1;
10820  }
10821  }
10822 
10823  /*
10824  if geometry is a point, a linestring or set of either and bounds not set,
10825  increase extent by a pixel to avoid missing points on border
10826 
10827  a whole pixel is used instead of half-pixel due to backward
10828  compatibility with GDAL 1.6, 1.7 and 1.8. 1.9+ works fine with half-pixel.
10829  */
10830  wkbtype = wkbFlatten(OGR_G_GetGeometryType(src_geom));
10831  if ((
10832  (wkbtype == wkbPoint) ||
10833  (wkbtype == wkbMultiPoint) ||
10834  (wkbtype == wkbLineString) ||
10835  (wkbtype == wkbMultiLineString)
10836  ) &&
10837  _dim[0] == 0 &&
10838  _dim[1] == 0
10839  ) {
10840  int result;
10841  LWPOLY *epoly = NULL;
10842  LWGEOM *lwgeom = NULL;
10843  GEOSGeometry *egeom = NULL;
10844  GEOSGeometry *geom = NULL;
10845 
10846  RASTER_DEBUG(3, "Testing geometry is properly contained by extent");
10847 
10848  /*
10849  see if geometry is properly contained by extent
10850  all parts of geometry lies within extent
10851  */
10852 
10853  /* initialize GEOS */
10854  initGEOS(lwnotice, lwgeom_geos_error);
10855 
10856  /* convert envelope to geometry */
10857  RASTER_DEBUG(4, "Converting envelope to geometry");
10858  epoly = rt_util_envelope_to_lwpoly(extent);
10859  if (epoly == NULL) {
10860  rterror("rt_raster_gdal_rasterize: Could not create envelope's geometry to test if geometry is properly contained by extent");
10861 
10863 
10864  OGR_G_DestroyGeometry(src_geom);
10865  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
10866  /* OGRCleanupAll(); */
10867 
10868  return NULL;
10869  }
10870 
10871  egeom = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(epoly));
10872  lwpoly_free(epoly);
10873 
10874  /* convert WKB to geometry */
10875  RASTER_DEBUG(4, "Converting WKB to geometry");
10876  lwgeom = lwgeom_from_wkb(wkb, wkb_len, LW_PARSER_CHECK_NONE);
10877  geom = (GEOSGeometry *) LWGEOM2GEOS(lwgeom);
10878  lwgeom_free(lwgeom);
10879 
10880  result = GEOSRelatePattern(egeom, geom, "T**FF*FF*");
10881  GEOSGeom_destroy(geom);
10882  GEOSGeom_destroy(egeom);
10883 
10884  if (result == 2) {
10885  rterror("rt_raster_gdal_rasterize: Could not test if geometry is properly contained by extent for geometry within extent");
10886 
10888 
10889  OGR_G_DestroyGeometry(src_geom);
10890  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
10891  /* OGRCleanupAll(); */
10892 
10893  return NULL;
10894  }
10895 
10896  /* geometry NOT properly contained by extent */
10897  if (!result) {
10898 
10899 #if POSTGIS_GDAL_VERSION > 18
10900 
10901  /* check alignment flag: grid_xw */
10902  if (
10903  (NULL == ul_xw && NULL == ul_yw) &&
10904  (NULL != grid_xw && NULL != grid_xw) &&
10905  FLT_NEQ(*grid_xw, extent.MinX)
10906  ) {
10907  // do nothing
10908  RASTER_DEBUG(3, "Skipping extent adjustment on X-axis due to upcoming alignment");
10909  }
10910  else {
10911  RASTER_DEBUG(3, "Adjusting extent for GDAL > 1.8 by half the scale on X-axis");
10912  extent.MinX -= (_scale[0] / 2.);
10913  extent.MaxX += (_scale[0] / 2.);
10914  }
10915 
10916  /* check alignment flag: grid_yw */
10917  if (
10918  (NULL == ul_xw && NULL == ul_yw) &&
10919  (NULL != grid_xw && NULL != grid_xw) &&
10920  FLT_NEQ(*grid_yw, extent.MaxY)
10921  ) {
10922  // do nothing
10923  RASTER_DEBUG(3, "Skipping extent adjustment on Y-axis due to upcoming alignment");
10924  }
10925  else {
10926  RASTER_DEBUG(3, "Adjusting extent for GDAL > 1.8 by half the scale on Y-axis");
10927  extent.MinY -= (_scale[1] / 2.);
10928  extent.MaxY += (_scale[1] / 2.);
10929  }
10930 
10931 #else
10932 
10933  /* check alignment flag: grid_xw */
10934  if (
10935  (NULL == ul_xw && NULL == ul_yw) &&
10936  (NULL != grid_xw && NULL != grid_xw) &&
10937  FLT_NEQ(*grid_xw, extent.MinX)
10938  ) {
10939  // do nothing
10940  RASTER_DEBUG(3, "Skipping extent adjustment on X-axis due to upcoming alignment");
10941  }
10942  else {
10943  RASTER_DEBUG(3, "Adjusting extent for GDAL <= 1.8 by the scale on X-axis");
10944  extent.MinX -= _scale[0];
10945  extent.MaxX += _scale[0];
10946  }
10947 
10948 
10949  /* check alignment flag: grid_yw */
10950  if (
10951  (NULL == ul_xw && NULL == ul_yw) &&
10952  (NULL != grid_xw && NULL != grid_xw) &&
10953  FLT_NEQ(*grid_yw, extent.MaxY)
10954  ) {
10955  // do nothing
10956  RASTER_DEBUG(3, "Skipping extent adjustment on Y-axis due to upcoming alignment");
10957  }
10958  else {
10959  RASTER_DEBUG(3, "Adjusting extent for GDAL <= 1.8 by the scale on Y-axis");
10960  extent.MinY -= _scale[1];
10961  extent.MaxY += _scale[1];
10962  }
10963 
10964 #endif
10965 
10966  }
10967 
10968  RASTER_DEBUGF(3, "Adjusted extent: %f, %f, %f, %f",
10969  extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
10970 
10971  extent.UpperLeftX = extent.MinX;
10972  extent.UpperLeftY = extent.MaxY;
10973  }
10974 
10975  /* reprocess extent if skewed */
10976  if (
10977  FLT_NEQ(_skew[0], 0) ||
10978  FLT_NEQ(_skew[1], 0)
10979  ) {
10980  rt_raster skewedrast;
10981 
10982  RASTER_DEBUG(3, "Computing skewed extent's envelope");
10983 
10984  skewedrast = rt_raster_compute_skewed_raster(
10985  extent,
10986  _skew,
10987  _scale,
10988  0.01
10989  );
10990  if (skewedrast == NULL) {
10991  rterror("rt_raster_gdal_rasterize: Could not compute skewed raster");
10992 
10994 
10995  OGR_G_DestroyGeometry(src_geom);
10996  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
10997  /* OGRCleanupAll(); */
10998 
10999  return NULL;
11000  }
11001 
11002  _dim[0] = skewedrast->width;
11003  _dim[1] = skewedrast->height;
11004 
11005  extent.UpperLeftX = skewedrast->ipX;
11006  extent.UpperLeftY = skewedrast->ipY;
11007 
11008  rt_raster_destroy(skewedrast);
11009  }
11010 
11011  /* raster dimensions */
11012  if (!_dim[0])
11013  _dim[0] = (int) fmax((fabs(extent.MaxX - extent.MinX) + (_scale[0] / 2.)) / _scale[0], 1);
11014  if (!_dim[1])
11015  _dim[1] = (int) fmax((fabs(extent.MaxY - extent.MinY) + (_scale[1] / 2.)) / _scale[1], 1);
11016 
11017  /* temporary raster */
11018  rast = rt_raster_new(_dim[0], _dim[1]);
11019  if (rast == NULL) {
11020  rterror("rt_raster_gdal_rasterize: Out of memory allocating temporary raster");
11021 
11023 
11024  OGR_G_DestroyGeometry(src_geom);
11025  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11026  /* OGRCleanupAll(); */
11027 
11028  return NULL;
11029  }
11030 
11031  /* set raster's spatial attributes */
11032  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
11033  rt_raster_set_scale(rast, _scale[0], -1 * _scale[1]);
11034  rt_raster_set_skews(rast, _skew[0], _skew[1]);
11035 
11037  RASTER_DEBUGF(3, "Temp raster's geotransform: %f, %f, %f, %f, %f, %f",
11038  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
11039  RASTER_DEBUGF(3, "Temp raster's dimensions (width x height): %d x %d",
11040  _dim[0], _dim[1]);
11041 
11042  /* user-specified upper-left corner */
11043  if (
11044  NULL != ul_xw &&
11045  NULL != ul_yw
11046  ) {
11047  ul_user = 1;
11048 
11049  RASTER_DEBUGF(4, "Using user-specified upper-left corner: %f, %f", *ul_xw, *ul_yw);
11050 
11051  /* set upper-left corner */
11052  rt_raster_set_offsets(rast, *ul_xw, *ul_yw);
11053  extent.UpperLeftX = *ul_xw;
11054  extent.UpperLeftY = *ul_yw;
11055  }
11056  else if (
11057  ((NULL != ul_xw) && (NULL == ul_yw)) ||
11058  ((NULL == ul_xw) && (NULL != ul_yw))
11059  ) {
11060  rterror("rt_raster_gdal_rasterize: Both X and Y upper-left corner values must be provided");
11061 
11062  rt_raster_destroy(rast);
11064 
11065  OGR_G_DestroyGeometry(src_geom);
11066  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11067  /* OGRCleanupAll(); */
11068 
11069  return NULL;
11070  }
11071 
11072  /* alignment only considered if upper-left corner not provided */
11073  if (
11074  !ul_user && (
11075  (NULL != grid_xw) || (NULL != grid_yw)
11076  )
11077  ) {
11078 
11079  if (
11080  ((NULL != grid_xw) && (NULL == grid_yw)) ||
11081  ((NULL == grid_xw) && (NULL != grid_yw))
11082  ) {
11083  rterror("rt_raster_gdal_rasterize: Both X and Y alignment values must be provided");
11084 
11085  rt_raster_destroy(rast);
11087 
11088  OGR_G_DestroyGeometry(src_geom);
11089  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11090  /* OGRCleanupAll(); */
11091 
11092  return NULL;
11093  }
11094 
11095  RASTER_DEBUGF(4, "Aligning extent to user-specified grid: %f, %f", *grid_xw, *grid_yw);
11096 
11097  do {
11098  double _r[2] = {0};
11099  double _w[2] = {0};
11100 
11101  /* raster is already aligned */
11102  if (FLT_EQ(*grid_xw, extent.UpperLeftX) && FLT_EQ(*grid_yw, extent.UpperLeftY)) {
11103  RASTER_DEBUG(3, "Skipping raster alignment as it is already aligned to grid");
11104  break;
11105  }
11106 
11107  extent.UpperLeftX = rast->ipX;
11108  extent.UpperLeftY = rast->ipY;
11109  rt_raster_set_offsets(rast, *grid_xw, *grid_yw);
11110 
11111  /* process upper-left corner */
11113  rast,
11114  extent.UpperLeftX, extent.UpperLeftY,
11115  &(_r[0]), &(_r[1]),
11116  NULL
11117  ) != ES_NONE) {
11118  rterror("rt_raster_gdal_rasterize: Could not compute raster pixel for spatial coordinates");
11119 
11120  rt_raster_destroy(rast);
11122 
11123  OGR_G_DestroyGeometry(src_geom);
11124  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11125  /* OGRCleanupAll(); */
11126 
11127  return NULL;
11128  }
11129 
11131  rast,
11132  _r[0], _r[1],
11133  &(_w[0]), &(_w[1]),
11134  NULL
11135  ) != ES_NONE) {
11136  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
11137 
11138  rt_raster_destroy(rast);
11140 
11141  OGR_G_DestroyGeometry(src_geom);
11142  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11143  /* OGRCleanupAll(); */
11144 
11145  return NULL;
11146  }
11147 
11148  /* shift occurred */
11149  if (FLT_NEQ(_w[0], extent.UpperLeftX)) {
11150  if (NULL == width)
11151  rast->width++;
11152  else if (NULL == scale_x) {
11153  double _c[2] = {0};
11154 
11155  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
11156 
11157  /* get upper-right corner */
11159  rast,
11160  rast->width, 0,
11161  &(_c[0]), &(_c[1]),
11162  NULL
11163  ) != ES_NONE) {
11164  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
11165 
11166  rt_raster_destroy(rast);
11168 
11169  OGR_G_DestroyGeometry(src_geom);
11170  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11171  /* OGRCleanupAll(); */
11172 
11173  return NULL;
11174  }
11175 
11176  rast->scaleX = fabs((_c[0] - _w[0]) / ((double) rast->width));
11177  }
11178  }
11179  if (FLT_NEQ(_w[1], extent.UpperLeftY)) {
11180  if (NULL == height)
11181  rast->height++;
11182  else if (NULL == scale_y) {
11183  double _c[2] = {0};
11184 
11185  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
11186 
11187  /* get upper-right corner */
11189  rast,
11190  0, rast->height,
11191  &(_c[0]), &(_c[1]),
11192  NULL
11193  ) != ES_NONE) {
11194  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
11195 
11196  rt_raster_destroy(rast);
11198 
11199  OGR_G_DestroyGeometry(src_geom);
11200  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11201  /* OGRCleanupAll(); */
11202 
11203  return NULL;
11204  }
11205 
11206  rast->scaleY = -1 * fabs((_c[1] - _w[1]) / ((double) rast->height));
11207  }
11208  }
11209 
11210  rt_raster_set_offsets(rast, _w[0], _w[1]);
11211  }
11212  while (0);
11213  }
11214 
11215  /*
11216  after this point, rt_envelope extent is no longer used
11217  */
11218 
11219  /* get key attributes from rast */
11220  _dim[0] = rast->width;
11221  _dim[1] = rast->height;
11223 
11224  /* scale-x is negative or scale-y is positive */
11225  if ((
11226  (NULL != scale_x) && (*scale_x < 0.)
11227  ) || (
11228  (NULL != scale_y) && (*scale_y > 0)
11229  )) {
11230  double _w[2] = {0};
11231 
11232  /* negative scale-x */
11233  if (
11234  (NULL != scale_x) &&
11235  (*scale_x < 0.)
11236  ) {
11237  RASTER_DEBUG(3, "Processing negative scale-x");
11238 
11240  rast,
11241  _dim[0], 0,
11242  &(_w[0]), &(_w[1]),
11243  NULL
11244  ) != ES_NONE) {
11245  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
11246 
11247  rt_raster_destroy(rast);
11249 
11250  OGR_G_DestroyGeometry(src_geom);
11251  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11252  /* OGRCleanupAll(); */
11253 
11254  return NULL;
11255  }
11256 
11257  _gt[0] = _w[0];
11258  _gt[1] = *scale_x;
11259 
11260  /* check for skew */
11261  if (NULL != skew_x && FLT_NEQ(*skew_x, 0))
11262  _gt[2] = *skew_x;
11263  }
11264  /* positive scale-y */
11265  if (
11266  (NULL != scale_y) &&
11267  (*scale_y > 0)
11268  ) {
11269  RASTER_DEBUG(3, "Processing positive scale-y");
11270 
11272  rast,
11273  0, _dim[1],
11274  &(_w[0]), &(_w[1]),
11275  NULL
11276  ) != ES_NONE) {
11277  rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
11278 
11279  rt_raster_destroy(rast);
11281 
11282  OGR_G_DestroyGeometry(src_geom);
11283  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11284  /* OGRCleanupAll(); */
11285 
11286  return NULL;
11287  }
11288 
11289  _gt[3] = _w[1];
11290  _gt[5] = *scale_y;
11291 
11292  /* check for skew */
11293  if (NULL != skew_y && FLT_NEQ(*skew_y, 0))
11294  _gt[4] = *skew_y;
11295  }
11296  }
11297 
11298  rt_raster_destroy(rast);
11299  rast = NULL;
11300 
11301  RASTER_DEBUGF(3, "Applied geotransform: %f, %f, %f, %f, %f, %f",
11302  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
11303  RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d",
11304  _dim[0], _dim[1]);
11305 
11306  /* load GDAL mem */
11307  if (!rt_util_gdal_driver_registered("MEM")) {
11308  GDALRegister_MEM();
11309  unload_drv = 1;
11310  }
11311  _drv = GDALGetDriverByName("MEM");
11312  if (NULL == _drv) {
11313  rterror("rt_raster_gdal_rasterize: Could not load the MEM GDAL driver");
11314 
11316 
11317  OGR_G_DestroyGeometry(src_geom);
11318  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11319  /* OGRCleanupAll(); */
11320 
11321  return NULL;
11322  }
11323 
11324  if (unload_drv)
11325  GDALDeregisterDriver(_drv);
11326 
11327  _ds = GDALCreate(_drv, "", _dim[0], _dim[1], 0, GDT_Byte, NULL);
11328  if (NULL == _ds) {
11329  rterror("rt_raster_gdal_rasterize: Could not create a GDALDataset to rasterize the geometry into");
11330 
11332 
11333  OGR_G_DestroyGeometry(src_geom);
11334  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11335  /* OGRCleanupAll(); */
11336 
11337  if (unload_drv) GDALDestroyDriver(_drv);
11338 
11339  return NULL;
11340  }
11341 
11342  /* set geotransform */
11343  cplerr = GDALSetGeoTransform(_ds, _gt);
11344  if (cplerr != CE_None) {
11345  rterror("rt_raster_gdal_rasterize: Could not set geotransform on GDALDataset");
11346 
11348 
11349  OGR_G_DestroyGeometry(src_geom);
11350  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11351  /* OGRCleanupAll(); */
11352 
11353  GDALClose(_ds);
11354  if (unload_drv) GDALDestroyDriver(_drv);
11355 
11356  return NULL;
11357  }
11358 
11359  /* set SRS */
11360  if (NULL != src_sr) {
11361  char *_srs = NULL;
11362  OSRExportToWkt(src_sr, &_srs);
11363 
11364  cplerr = GDALSetProjection(_ds, _srs);
11365  CPLFree(_srs);
11366  if (cplerr != CE_None) {
11367  rterror("rt_raster_gdal_rasterize: Could not set projection on GDALDataset");
11368 
11370 
11371  OGR_G_DestroyGeometry(src_geom);
11372  OSRDestroySpatialReference(src_sr);
11373  /* OGRCleanupAll(); */
11374 
11375  GDALClose(_ds);
11376  if (unload_drv) GDALDestroyDriver(_drv);
11377 
11378  return NULL;
11379  }
11380  }
11381 
11382  /* set bands */
11383  for (i = 0; i < arg->numbands; i++) {
11384  err = 0;
11385 
11386  do {
11387  /* add band */
11388  cplerr = GDALAddBand(_ds, rt_util_pixtype_to_gdal_datatype(arg->pixtype[i]), NULL);
11389  if (cplerr != CE_None) {
11390  rterror("rt_raster_gdal_rasterize: Could not add band to GDALDataset");
11391  err = 1;
11392  break;
11393  }
11394 
11395  _band = GDALGetRasterBand(_ds, i + 1);
11396  if (NULL == _band) {
11397  rterror("rt_raster_gdal_rasterize: Could not get band %d from GDALDataset", i + 1);
11398  err = 1;
11399  break;
11400  }
11401 
11402  /* nodata value */
11403  if (arg->hasnodata[i]) {
11404  RASTER_DEBUGF(4, "Setting NODATA value of band %d to %f", i, arg->nodata[i]);
11405  cplerr = GDALSetRasterNoDataValue(_band, arg->nodata[i]);
11406  if (cplerr != CE_None) {
11407  rterror("rt_raster_gdal_rasterize: Could not set nodata value");
11408  err = 1;
11409  break;
11410  }
11411  RASTER_DEBUGF(4, "NODATA value set to %f", GDALGetRasterNoDataValue(_band, NULL));
11412  }
11413 
11414  /* initial value */
11415  cplerr = GDALFillRaster(_band, arg->init[i], 0);
11416  if (cplerr != CE_None) {
11417  rterror("rt_raster_gdal_rasterize: Could not set initial value");
11418  err = 1;
11419  break;
11420  }
11421  }
11422  while (0);
11423 
11424  if (err) {
11426 
11427  OGR_G_DestroyGeometry(src_geom);
11428  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11429  /* OGRCleanupAll(); */
11430 
11431  GDALClose(_ds);
11432  if (unload_drv) GDALDestroyDriver(_drv);
11433 
11434  return NULL;
11435  }
11436  }
11437 
11438  arg->bandlist = (int *) rtalloc(sizeof(int) * arg->numbands);
11439  for (i = 0; i < arg->numbands; i++) arg->bandlist[i] = i + 1;
11440 
11441 
11442  /* burn geometry */
11443  cplerr = GDALRasterizeGeometries(
11444  _ds,
11445  arg->numbands, arg->bandlist,
11446  1, &src_geom,
11447  NULL, NULL,
11448  arg->value,
11449  options,
11450  NULL, NULL
11451  );
11452  if (cplerr != CE_None) {
11453  rterror("rt_raster_gdal_rasterize: Could not rasterize geometry");
11454 
11456 
11457  OGR_G_DestroyGeometry(src_geom);
11458  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11459  /* OGRCleanupAll(); */
11460 
11461  GDALClose(_ds);
11462  if (unload_drv) GDALDestroyDriver(_drv);
11463 
11464  return NULL;
11465  }
11466 
11467  /* convert gdal dataset to raster */
11468  GDALFlushCache(_ds);
11469  RASTER_DEBUG(3, "Converting GDAL dataset to raster");
11470  rast = rt_raster_from_gdal_dataset(_ds);
11471 
11472  OGR_G_DestroyGeometry(src_geom);
11473  if (src_sr != NULL) OSRDestroySpatialReference(src_sr);
11474  /* OGRCleanupAll(); */
11475 
11476  GDALClose(_ds);
11477  if (unload_drv) GDALDestroyDriver(_drv);
11478 
11479  if (NULL == rast) {
11480  rterror("rt_raster_gdal_rasterize: Could not rasterize geometry");
11481  return NULL;
11482  }
11483 
11484  /* width, height */
11485  _width = rt_raster_get_width(rast);
11486  _height = rt_raster_get_height(rast);
11487 
11488  /* check each band for pixtype */
11489  for (i = 0; i < arg->numbands; i++) {
11490  uint8_t *data = NULL;
11491  rt_band band = NULL;
11492  rt_band oldband = NULL;
11493 
11494  double val = 0;
11495  int nodata = 0;
11496  int hasnodata = 0;
11497  double nodataval = 0;
11498  int x = 0;
11499  int y = 0;
11500 
11501  oldband = rt_raster_get_band(rast, i);
11502  if (oldband == NULL) {
11503  rterror("rt_raster_gdal_rasterize: Could not get band %d of output raster", i);
11505  rt_raster_destroy(rast);
11506  return NULL;
11507  }
11508 
11509  /* band is of user-specified type */
11510  if (rt_band_get_pixtype(oldband) == arg->pixtype[i])
11511  continue;
11512 
11513  /* hasnodata, nodataval */
11514  hasnodata = rt_band_get_hasnodata_flag(oldband);
11515  if (hasnodata)
11516  rt_band_get_nodata(oldband, &nodataval);
11517 
11518  /* allocate data */
11519  data = rtalloc(rt_pixtype_size(arg->pixtype[i]) * _width * _height);
11520  if (data == NULL) {
11521  rterror("rt_raster_gdal_rasterize: Could not allocate memory for band data");
11523  rt_raster_destroy(rast);
11524  return NULL;
11525  }
11526  memset(data, 0, rt_pixtype_size(arg->pixtype[i]) * _width * _height);
11527 
11528  /* create new band of correct type */
11529  band = rt_band_new_inline(
11530  _width, _height,
11531  arg->pixtype[i],
11532  hasnodata, nodataval,
11533  data
11534  );
11535  if (band == NULL) {
11536  rterror("rt_raster_gdal_rasterize: Could not create band");
11537  rtdealloc(data);
11539  rt_raster_destroy(rast);
11540  return NULL;
11541  }
11542 
11543  /* give ownership of data to band */
11544  rt_band_set_ownsdata_flag(band, 1);
11545 
11546  /* copy pixel by pixel */
11547  for (x = 0; x < _width; x++) {
11548  for (y = 0; y < _height; y++) {
11549  err = rt_band_get_pixel(oldband, x, y, &val, &nodata);
11550  if (err != ES_NONE) {
11551  rterror("rt_raster_gdal_rasterize: Could not get pixel value");
11553  rt_raster_destroy(rast);
11554  rt_band_destroy(band);
11555  return NULL;
11556  }
11557 
11558  if (nodata)
11559  val = nodataval;
11560 
11561  err = rt_band_set_pixel(band, x, y, val, NULL);
11562  if (err != ES_NONE) {
11563  rterror("rt_raster_gdal_rasterize: Could not set pixel value");
11565  rt_raster_destroy(rast);
11566  rt_band_destroy(band);
11567  return NULL;
11568  }
11569  }
11570  }
11571 
11572  /* replace band */
11573  oldband = rt_raster_replace_band(rast, band, i);
11574  if (oldband == NULL) {
11575  rterror("rt_raster_gdal_rasterize: Could not replace band %d of output raster", i);
11577  rt_raster_destroy(rast);
11578  rt_band_destroy(band);
11579  return NULL;
11580  }
11581 
11582  /* free oldband */
11583  rt_band_destroy(oldband);
11584  }
11585 
11587 
11588  RASTER_DEBUG(3, "done");
11589 
11590  return rast;
11591 }
rt_raster rt_raster_from_gdal_dataset(GDALDatasetH ds)
Return a raster from a GDAL dataset.
Definition: rt_api.c:9339
#define FLT_NEQ(x, y)
Definition: rt_api.h:2158
double MinY
Definition: rt_api.h:154
void rt_util_from_ogr_envelope(OGREnvelope env, rt_envelope *ext)
Definition: rt_api.c:514
void rtdealloc(void *mem)
Definition: rt_api.c:882
double UpperLeftX
Definition: rt_api.h:157
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_api.c:5387
rt_raster rt_raster_compute_skewed_raster(rt_envelope extent, double *skew, double *scale, double tolerance)
Definition: rt_api.c:6764
tuple data
Definition: ovdump.py:103
double MaxY
Definition: rt_api.h:155
double MaxX
Definition: rt_api.h:153
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1006
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_api.c:1900
rt_band rt_raster_replace_band(rt_raster raster, rt_band band, int index)
Replace band at provided index with new band.
Definition: rt_api.c:8702
tuple band
Definition: ovdump.py:57
tuple rast
Definition: rtpixdump.py:62
double ipY
Definition: rt_api.h:2221
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom)
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_api.c:1466
char ** result
Definition: liblwgeom.h:218
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:239
uint16_t height
Definition: rt_api.h:2227
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
Definition: rt_api.c:5473
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition: rt_api.c:1097
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_api.c:6054
void rtinfo(const char *fmt,...)
Definition: rt_api.c:907
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_api.c:6005
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:1706
#define RASTER_DEBUG(level, msg)
Definition: rt_api.h:281
int rt_util_gdal_driver_registered(const char *drv)
Definition: rt_api.c:461
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:54
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_api.c:5442
rt_pixtype
Definition: rt_api.h:172
#define RASTER_DEBUGF(level, msg,...)
Definition: rt_api.h:285
void lwgeom_geos_error(const char *fmt,...)
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:79
void rt_band_set_ownsdata_flag(rt_band band, int flag)
Definition: rt_api.c:1936
uint16_t width
Definition: rt_api.h:2226
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_api.c:3058
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_api.c:2549
rt_pixtype * pixtype
Definition: rt_api.c:10541
static _rti_rasterize_arg _rti_rasterize_arg_init()
Definition: rt_api.c:10550
#define FLT_EQ(x, y)
Definition: rt_api.h:2159
tuple x
Definition: pixval.py:53
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_api.c:1650
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_api.c:5434
void * rtalloc(size_t size)
Raster core memory management functions.
Definition: rt_api.c:867
static void _rti_rasterize_arg_destroy(_rti_rasterize_arg arg)
Definition: rt_api.c:10573
void rterror(const char *fmt,...)
Raster core error and info handlers.
Definition: rt_api.c:895
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_api.c:5686
double MinX
Definition: rt_api.h:152
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_api.c:6105
double ipX
Definition: rt_api.h:2220
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_api.c:5353
double scaleX
Definition: rt_api.h:2218
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:728
LWPOLY * rt_util_envelope_to_lwpoly(rt_envelope env)
Definition: rt_api.c:543
GDALDataType rt_util_pixtype_to_gdal_datatype(rt_pixtype pt)
Convert rt_pixtype to GDALDataType.
Definition: rt_api.c:229
double UpperLeftY
Definition: rt_api.h:158
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_api.c:2002
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_api.c:5504
tuple y
Definition: pixval.py:54
double scaleY
Definition: rt_api.h:2219
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_api.c:5426
uint8_t * hasnodata
Definition: rt_api.c:10544
rt_errorstate rt_band_set_pixel(rt_band band, int x, int y, double val, int *converted)
Set single pixel's value.
Definition: rt_api.c:2302

Here is the call graph for this function:

Here is the caller graph for this function: