GDALTableFunctions Namespace Reference

Functions
template<typename TV , typename TG >
int32_t	tf_raster_contour_impl (TableFunctionManager &mgr, const int32_t raster_width, const int32_t raster_height, const double affine_transform[6], const TV *values, const TV contour_interval, const TV contour_offset, Column< TG > &contour_features, Column< TV > &contour_values)
template<typename TLL , typename TV , typename TG >
int32_t	tf_raster_contour_rasterize_impl (TableFunctionManager &mgr, const Column< TLL > &lon, const Column< TLL > &lat, const Column< TV > &values, const TextEncodingNone &agg_type, const float bin_dim_meters, const int32_t neighborhood_fill_radius, const bool fill_only_nulls, const TextEncodingNone &fill_agg_type, const bool flip_latitude, const TV contour_interval, const TV contour_offset, Column< TG > &contour_features, Column< TV > &contour_values)
template<typename TLL , typename TV , typename TG >
int32_t	tf_raster_contour_direct_impl (TableFunctionManager &mgr, const Column< TLL > &lon, const Column< TLL > &lat, const Column< TV > &values, const int32_t raster_width, const int32_t raster_height, const bool flip_latitude, const TV contour_interval, const TV contour_offset, Column< TG > &contour_features, Column< TV > &contour_values)

Function Documentation

template<typename TLL , typename TV , typename TG >

int32_t GDALTableFunctions::tf_raster_contour_direct_impl	(	TableFunctionManager &	mgr,
		const Column< TLL > &	lon,
		const Column< TLL > &	lat,
		const Column< TV > &	values,
		const int32_t	raster_width,
		const int32_t	raster_height,
		const bool	flip_latitude,
		const TV	contour_interval,
		const TV	contour_offset,
		Column< TG > &	contour_features,
		Column< TV > &	contour_values
	)

Definition at line 485 of file GDALTableFunctions.hpp.

References CHECK_GE, CHECK_GT, Column< T >::getPtr(), and Column< T >::size().

                                                                  {
  // enforce parser validation
  CHECK_GE(raster_width, 1u);
  CHECK_GE(raster_height, 1u);
  CHECK_GT(contour_interval, static_cast<TV>(0));

  // expected pixel counts
  auto const num_pixels =
      static_cast<int64_t>(raster_width) * static_cast<int64_t>(raster_height);
  if (lon.size() != num_pixels || lat.size() != num_pixels ||
      values.size() != num_pixels) {
    return mgr.ERROR_MESSAGE("Raster lon/lat/values size mismatch");
  }

  // find lon/lat range
  double lon_min = double(*std::min_element(lon.getPtr(), lon.getPtr() + num_pixels));
  double lon_max = double(*std::max_element(lon.getPtr(), lon.getPtr() + num_pixels));
  double lat_min = double(*std::min_element(lat.getPtr(), lat.getPtr() + num_pixels));
  double lat_max = double(*std::max_element(lat.getPtr(), lat.getPtr() + num_pixels));

  // build affine transform matrix
  // @TODO(se) verify the -1s
  std::array<double, 6> affine_transform;
  affine_transform[0] = lon_min;
  affine_transform[1] = (lon_max - lon_min) / double(raster_width - 1);
  affine_transform[2] = 0.0;
  affine_transform[4] = 0.0;
  if (flip_latitude) {
    affine_transform[3] = lat_max;
    affine_transform[5] = (lat_min - lat_max) / double(raster_height - 1);
  } else {
    affine_transform[3] = lat_min;
    affine_transform[5] = (lat_max - lat_min) / double(raster_height - 1);
  }

  // do it
  return tf_raster_contour_impl<TV, TG>(mgr,
                                        raster_width,
                                        raster_height,
                                        affine_transform.data(),
                                        values.getPtr(),
                                        contour_interval,
                                        contour_offset,
                                        contour_features,
                                        contour_values);
}

Here is the call graph for this function:

template<typename TV , typename TG >

int32_t GDALTableFunctions::tf_raster_contour_impl	(	TableFunctionManager &	mgr,
		const int32_t	raster_width,
		const int32_t	raster_height,
		const double	affine_transform[6],
		const TV *	values,
		const TV	contour_interval,
		const TV	contour_offset,
		Column< TG > &	contour_features,
		Column< TV > &	contour_values
	)

Definition at line 39 of file GDALTableFunctions.hpp.

References CHECK, CHECK_GE, Geospatial::GDAL::init(), run_benchmark_import::result, TableFunctionManager::set_output_array_values_total_number(), TableFunctionManager::set_output_row_size(), toString(), and VLOG.

                                                           {
  // supported types
  static_assert(std::is_same_v<float, TV> || std::is_same_v<double, TV>);
  static_assert(std::is_same_v<GeoLineString, TG> || std::is_same_v<GeoPolygon, TG>);

  // value type
  std::string value_type_str;
  int pixel_offset{};
  double null_value;
  if constexpr (std::is_same_v<float, TV>) {
    value_type_str = "Float32";
    pixel_offset = sizeof(float);
    null_value = static_cast<double>(std::numeric_limits<float>::min());
  } else if constexpr (std::is_same_v<double, TV>) {
    value_type_str = "Float64";
    pixel_offset = sizeof(double);
    null_value = std::numeric_limits<double>::min();
  }
  CHECK(value_type_str.length());

  auto const num_values = raster_width * raster_height;
  VLOG(2) << "tf_raster_contour: Final raster data has " << num_values << " values, min "
          << *std::min_element(values, values + num_values) << ", max "
          << *std::max_element(values, values + num_values);

  // geo type
  constexpr bool is_polygons = std::is_same_v<GeoPolygon, TG>;

  // construct in-memory-raster string
  std::stringstream mem_ss;
  mem_ss << "MEM:::DATAPOINTER=" << std::hex << values << std::dec
         << ",PIXELS=" << raster_width << ",LINES=" << raster_height
         << ",BANDS=1,DATATYPE=" << value_type_str << ",PIXELOFFSET=" << pixel_offset
         << ",GEOTRANSFORM=" << affine_transform[0] << "/" << affine_transform[1] << "/"
         << affine_transform[2] << "/" << affine_transform[3] << "/"
         << affine_transform[4] << "/" << affine_transform[5];
  std::string mem_str = mem_ss.str();

  // lazy init GDAL
  Geospatial::GDAL::init();

  // things to clean up
  char** options = nullptr;
  GDALDatasetH raster_ds = nullptr;
  GDALDataset* vector_ds = nullptr;

  // auto clean-up on any exit
  ScopeGuard cleanup = [options, raster_ds, vector_ds]() {
    if (options) {
      CSLDestroy(options);
    }
    if (raster_ds) {
      GDALClose(raster_ds);
    }
    if (vector_ds) {
      GDALClose(vector_ds);
    }
  };

  // input dataset
  raster_ds = GDALOpen(mem_str.c_str(), GA_ReadOnly);
  CHECK(raster_ds);

  // input band
  auto raster_band = GDALGetRasterBand(raster_ds, 1);
  CHECK(raster_band);

  // output driver
  auto* memory_driver = GetGDALDriverManager()->GetDriverByName("Memory");
  CHECK(memory_driver);

  // output dataset
  vector_ds = memory_driver->Create("contours", 0, 0, 0, GDT_Unknown, NULL);
  CHECK(vector_ds);

  // output spatial reference
  OGRSpatialReference spatial_reference;
  spatial_reference.importFromEPSG(4326);
  spatial_reference.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

  // output layer
  auto* vector_layer =
      vector_ds->CreateLayer("lines", &spatial_reference, wkbLineString, NULL);
  CHECK(vector_layer);

  // contour values fields
  int contour_min_field_index{-1};
  int contour_max_field_index{-1};
  int contour_val_field_index{-1};
  if (is_polygons) {
    OGRFieldDefn contour_min_field_defn("contour_min", OFTReal);
    vector_layer->CreateField(&contour_min_field_defn);
    contour_min_field_index = vector_layer->FindFieldIndex("contour_min", TRUE);
    CHECK_GE(contour_min_field_index, 0);
    OGRFieldDefn contour_max_field_defn("contour_max", OFTReal);
    vector_layer->CreateField(&contour_max_field_defn);
    contour_max_field_index = vector_layer->FindFieldIndex("contour_max", TRUE);
    CHECK_GE(contour_max_field_index, 0);
  } else {
    OGRFieldDefn contour_val_field_defn("contour_val", OFTReal);
    vector_layer->CreateField(&contour_val_field_defn);
    contour_val_field_index = vector_layer->FindFieldIndex("contour_val", TRUE);
    CHECK_GE(contour_val_field_index, 0);
  }

  // options
  if constexpr (is_polygons) {
    options = CSLAppendPrintf(options, "ELEV_FIELD_MIN=%d", contour_min_field_index);
    options = CSLAppendPrintf(options, "ELEV_FIELD_MAX=%d", contour_max_field_index);
  } else {
    options = CSLAppendPrintf(options, "ELEV_FIELD=%d", contour_val_field_index);
  }
  options =
      CSLAppendPrintf(options, "LEVEL_INTERVAL=%f", static_cast<float>(contour_interval));
  options = CSLAppendPrintf(options, "LEVEL_BASE=%f", static_cast<float>(contour_offset));
  options = CSLAppendPrintf(options, "NODATA=%.19g", null_value);
  if constexpr (is_polygons) {
    options = CSLAppendPrintf(options, "POLYGONIZE=YES");
  }

  // generate contour
  auto result =
      GDALContourGenerateEx(raster_band, vector_layer, options, nullptr, nullptr);
  if (result != CE_None) {
    return mgr.ERROR_MESSAGE("Contour generation failed");
  }

  // reset the output dataset
  vector_ds->ResetReading();

  // get feature count
  auto const num_features = static_cast<int32_t>(vector_layer->GetFeatureCount());

  VLOG(1) << "tf_raster_contour: GDAL generated " << num_features << " features";

  // did we get any features?
  if (num_features == 0) {
    mgr.set_output_array_values_total_number(0, 0);
    mgr.set_output_array_values_total_number(1, 0);
    mgr.set_output_row_size(0);
    return 0;
  }

  // first pass, accumulate total sizes

  int64_t total_num_points{};
  int32_t num_output_features{};

  vector_layer->ResetReading();

  // iterate features
  for (int32_t i = 0; i < num_features; i++) {
    // get feature
    Geospatial::GDAL::FeatureUqPtr feature(vector_layer->GetNextFeature());

    // get geometry
    auto const* geometry = feature->GetGeometryRef();
    CHECK(geometry);

    // get geometry type
    auto const geometry_wkb_type = wkbFlatten(geometry->getGeometryType());

    if constexpr (is_polygons) {
      // check type
      if (geometry_wkb_type != wkbMultiPolygon) {
        return mgr.ERROR_MESSAGE("Geometry WKB type is not MultiPolygon");
      }

      // it's a polygon
      auto const* multipolygon_geometry = dynamic_cast<const OGRMultiPolygon*>(geometry);
      CHECK(multipolygon_geometry);

      // count the points from all rings of all polygons
      auto const num_geometries = multipolygon_geometry->getNumGeometries();
      for (auto p = 0; p < num_geometries; p++) {
        // get this polygon
        auto const* polygon_geometry =
            dynamic_cast<const OGRPolygon*>(multipolygon_geometry->getGeometryRef(p));
        CHECK(polygon_geometry);

        // count points in exterior ring
        auto const* exterior_ring = polygon_geometry->getExteriorRing();
        CHECK(exterior_ring);
        total_num_points += exterior_ring->getNumPoints();

        // count points in interior rings
        for (auto r = 0; r < polygon_geometry->getNumInteriorRings(); r++) {
          auto const* interior_ring = polygon_geometry->getInteriorRing(r);
          CHECK(interior_ring);
          total_num_points += interior_ring->getNumPoints();
        }

        // each polygon is an output feature
        num_output_features++;
      }
    } else {
      // check type
      if (geometry_wkb_type != wkbLineString) {
        return mgr.ERROR_MESSAGE("Geometry WKB type is not Linestring");
      }

      // it's a linestring
      auto const* linestring_geometry = dynamic_cast<const OGRLineString*>(geometry);
      CHECK(linestring_geometry);

      // accumulate these points
      total_num_points += linestring_geometry->getNumPoints();

      // one output feature
      num_output_features++;
    }
  }

  VLOG(1) << "tf_raster_contour: Total points " << total_num_points;

  // size outputs
  mgr.set_output_array_values_total_number(0, total_num_points * 2);
  mgr.set_output_array_values_total_number(1, num_output_features);
  mgr.set_output_row_size(num_output_features);

  // second pass, build output features

  int64_t output_feature_index{};

  vector_layer->ResetReading();

  // iterate features
  for (int32_t i = 0; i < num_features; i++) {
    // get feature
    Geospatial::GDAL::FeatureUqPtr feature(vector_layer->GetNextFeature());

    // get geometry
    auto const* geometry = feature->GetGeometryRef();
    CHECK(geometry);

    if constexpr (is_polygons) {
      // it's a polygon
      auto const* multipolygon_geometry = dynamic_cast<const OGRMultiPolygon*>(geometry);
      CHECK(multipolygon_geometry);

      // process each polygon as an output feature
      auto const num_geometries = multipolygon_geometry->getNumGeometries();
      for (auto p = 0; p < num_geometries; p++) {
        // get this polygon
        auto const* polygon_geometry =
            dynamic_cast<const OGRPolygon*>(multipolygon_geometry->getGeometryRef(p));
        CHECK(polygon_geometry);

        // gather the points from all rings
        std::vector<std::vector<double>> coords;
        {
          std::vector<double> ring_coords;
          auto const* exterior_ring = polygon_geometry->getExteriorRing();
          CHECK(exterior_ring);
          for (int j = 0; j < exterior_ring->getNumPoints(); j++) {
            OGRPoint point;
            exterior_ring->getPoint(j, &point);
            ring_coords.push_back(point.getX());
            ring_coords.push_back(point.getY());
          }
          coords.push_back(ring_coords);
        }
        for (auto r = 0; r < polygon_geometry->getNumInteriorRings(); r++) {
          auto const* interior_ring = polygon_geometry->getInteriorRing(r);
          CHECK(interior_ring);
          std::vector<double> ring_coords;
          for (int j = 0; j < interior_ring->getNumPoints(); j++) {
            OGRPoint point;
            interior_ring->getPoint(j, &point);
            ring_coords.push_back(point.getX());
            ring_coords.push_back(point.getY());
          }
          coords.push_back(ring_coords);
        }

        // set output contour polygon
        auto status = contour_features[output_feature_index].fromCoords(coords);
        if (status != FlatBufferManager::Status::Success) {
          return mgr.ERROR_MESSAGE("fromCoords failed: " + ::toString(status));
        }

        // set output contour value
        // min value may be zero; if so, replace with (max - interval)
        // that will still be valid in the cases of max == interval or offset != 0
        auto const contour_min =
            static_cast<TV>(feature->GetFieldAsDouble(contour_min_field_index));
        auto const contour_max =
            static_cast<TV>(feature->GetFieldAsDouble(contour_max_field_index));
        contour_values[output_feature_index] = (contour_min == static_cast<TV>(0))
                                                   ? (contour_max - contour_interval)
                                                   : contour_min;

        // done
        output_feature_index++;
      }
    } else {
      // it's a linestring
      auto const* linestring_geometry = dynamic_cast<const OGRLineString*>(geometry);
      CHECK(linestring_geometry);

      // unpack linestring
      std::vector<double> coords;
      for (int j = 0; j < linestring_geometry->getNumPoints(); j++) {
        OGRPoint point;
        linestring_geometry->getPoint(j, &point);
        coords.push_back(point.getX());
        coords.push_back(point.getY());
      }

      // set output contour linestring
      contour_features[output_feature_index].fromCoords(coords);

      // set output contour value
      contour_values[output_feature_index] =
          static_cast<TV>(feature->GetFieldAsDouble(contour_val_field_index));

      // done
      output_feature_index++;
    }
  }

  VLOG(1) << "tf_raster_contour: Output " << num_features << " features";

  // done
  return num_output_features;
}

Here is the call graph for this function:

template<typename TLL , typename TV , typename TG >

int32_t GDALTableFunctions::tf_raster_contour_rasterize_impl	(	TableFunctionManager &	mgr,
		const Column< TLL > &	lon,
		const Column< TLL > &	lat,
		const Column< TV > &	values,
		const TextEncodingNone &	agg_type,
		const float	bin_dim_meters,
		const int32_t	neighborhood_fill_radius,
		const bool	fill_only_nulls,
		const TextEncodingNone &	fill_agg_type,
		const bool	flip_latitude,
		const TV	contour_interval,
		const TV	contour_offset,
		Column< TG > &	contour_features,
		Column< TV > &	contour_values
	)

Definition at line 379 of file GDALTableFunctions.hpp.

References CHECK, CHECK_GE, CHECK_GT, f(), get_raster_agg_type(), Column< T >::getPtr(), TextEncodingNone::getString(), INVALID, TableFunctionManager::set_output_array_values_total_number(), TableFunctionManager::set_output_row_size(), Column< T >::size(), and VLOG.

                                                                     {
  // enforce parser validation
  CHECK_GT(bin_dim_meters, 0.0f);
  CHECK_GE(neighborhood_fill_radius, 0);
  CHECK_GT(contour_interval, static_cast<TV>(0));

  // validate the rest
  auto const raster_agg_type = get_raster_agg_type(agg_type, false);
  if (raster_agg_type == RasterAggType::INVALID) {
    auto const error_msg = "Invalid Raster Aggregate Type: " + agg_type.getString();
    return mgr.ERROR_MESSAGE(error_msg);
  }
  auto const raster_fill_agg_type = get_raster_agg_type(fill_agg_type, true);
  if (raster_fill_agg_type == RasterAggType::INVALID) {
    auto const error_msg =
        "Invalid Raster Fill Aggregate Type: " + fill_agg_type.getString();
    return mgr.ERROR_MESSAGE(error_msg);
  }

  VLOG(2) << "tf_raster_contour: Input raster data has " << values.size()
          << " values, min "
          << *std::min_element(values.getPtr(), values.getPtr() + values.size())
          << ", max "
          << *std::max_element(values.getPtr(), values.getPtr() + values.size());

  std::unique_ptr<GeoRaster<TLL, TV>> geo_raster;
  try {
    geo_raster = std::make_unique<GeoRaster<TLL, TV>>(
        lon, lat, values, raster_agg_type, bin_dim_meters, true, true);
    CHECK(geo_raster);
  } catch (std::runtime_error& e) {
    // the only exception this variant of GeoRaster constructor will throw
    // is when num_x_bins_ < 1 || num_y_bins_ < 1, which is handled below
    geo_raster = nullptr;
  }

  // GeoRaster constructor can still return zero size if no input
  if (!geo_raster || geo_raster->num_x_bins_ < 1 || geo_raster->num_y_bins_ < 1) {
    VLOG(1) << "tf_raster_contour: No input raster data. Cannot compute contours.";
    mgr.set_output_array_values_total_number(0, 0);
    mgr.set_output_array_values_total_number(1, 0);
    mgr.set_output_row_size(0);
    return 0;
  }

  if (neighborhood_fill_radius > 0) {
    geo_raster->fill_bins_from_neighbors(
        neighborhood_fill_radius, fill_only_nulls, raster_fill_agg_type);
  }

  auto const raster_width = static_cast<int32_t>(geo_raster->num_x_bins_);
  auto const raster_height = static_cast<int32_t>(geo_raster->num_y_bins_);
  auto const lon_min = static_cast<double>(geo_raster->x_min_);
  auto const lat_min = static_cast<double>(geo_raster->y_min_);
  auto const lon_bin_scale = static_cast<double>(geo_raster->x_scale_input_to_bin_);
  auto const lat_bin_scale = static_cast<double>(geo_raster->y_scale_input_to_bin_);

  // this should never happen, but make it an exception anyway
  if (lon_bin_scale <= 0.0 || lat_bin_scale <= 0.0) {
    return mgr.ERROR_MESSAGE("Invalid input raster scale. Cannot compute contours.");
  }

  // build affine transform matrix
  std::array<double, 6> affine_transform;
  affine_transform[0] = lon_min;
  affine_transform[1] = 1.0 / lon_bin_scale;
  affine_transform[2] = 0.0;
  affine_transform[4] = 0.0;
  if (flip_latitude) {
    affine_transform[3] = lat_min + (double(raster_height) / lat_bin_scale);
    affine_transform[5] = -1.0 / lat_bin_scale;
  } else {
    affine_transform[3] = lat_min;
    affine_transform[5] = 1.0 / lat_bin_scale;
  }

  // do it
  return tf_raster_contour_impl<TV, TG>(mgr,
                                        raster_width,
                                        raster_height,
                                        affine_transform.data(),
                                        geo_raster->z_.data(),
                                        contour_interval,
                                        contour_offset,
                                        contour_features,
                                        contour_values);
}

Here is the call graph for this function: