Kontinuation commented on code in PR #698: URL: https://github.com/apache/sedona-db/pull/698#discussion_r3001721353
########## c/sedona-gdal/src/raster/rasterize_affine.rs: ########## @@ -0,0 +1,362 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +//! Fast affine-transformer rasterize wrapper. +//! +//! GDALRasterizeGeometries() will internally call GDALCreateGenImgProjTransformer2() +//! if pfnTransformer is NULL, even in the common case where only a GeoTransform-based +//! affine conversion from georeferenced coords to pixel/line is needed. +//! +//! This module supplies a minimal GDALTransformerFunc that applies the dataset +//! GeoTransform (and its inverse), avoiding expensive transformer creation. + +use std::ffi::{c_int, c_void}; +use std::ptr; + +use crate::dataset::Dataset; +use crate::errors::{GdalError, Result}; +use crate::gdal_api::{call_gdal_api, GdalApi}; +use crate::gdal_dyn_bindgen::CE_None; +use crate::geo_transform::{GeoTransform, GeoTransformEx}; +use crate::vector::geometry::Geometry; + +#[repr(C)] +struct AffineTransformArg { + gt: GeoTransform, + inv_gt: GeoTransform, +} + +unsafe extern "C" fn affine_transformer( + p_transformer_arg: *mut c_void, + b_dst_to_src: c_int, + n_point_count: c_int, + x: *mut f64, + y: *mut f64, + _z: *mut f64, + pan_success: *mut c_int, +) -> c_int { + if p_transformer_arg.is_null() || x.is_null() || y.is_null() || pan_success.is_null() { + return 0; + } + if n_point_count < 0 { + return 0; + } + + // Treat transformer arg as immutable. + let arg = &*(p_transformer_arg as *const AffineTransformArg); + let affine = if b_dst_to_src == 0 { + &arg.inv_gt + } else { + &arg.gt + }; + + let n = n_point_count as usize; + for i in 0..n { + // SAFETY: x/y/pan_success are assumed to point to arrays of length n_point_count. + let xin = unsafe { *x.add(i) }; + let yin = unsafe { *y.add(i) }; + let (xout, yout) = affine.apply(xin, yin); + unsafe { + *x.add(i) = xout; + *y.add(i) = yout; + *pan_success.add(i) = 1; + } + } + + 1 +} + +/// Rasterize geometries using the dataset geotransform as the transformer. +/// Geometry coordinates must already be in the dataset georeferenced coordinate space. +pub fn rasterize_affine( + api: &'static GdalApi, + dataset: &Dataset, + bands: &[usize], + geometries: &[Geometry], + burn_values: &[f64], + all_touched: bool, +) -> Result<()> { + if bands.is_empty() { + return Err(GdalError::BadArgument( + "`bands` must not be empty".to_string(), + )); + } + if burn_values.len() != geometries.len() { + return Err(GdalError::BadArgument(format!( + "Burn values length ({}) must match geometries length ({})", + burn_values.len(), + geometries.len() + ))); + } + + let raster_count = dataset.raster_count(); + for band in bands { + let is_good = *band > 0 && *band <= raster_count; + if !is_good { + return Err(GdalError::BadArgument(format!( + "Band index {} is out of bounds", + *band + ))); + } + } + + let bands_i32: Vec<c_int> = bands.iter().map(|&band| band as c_int).collect(); + + let c_options = if all_touched { + [c"ALL_TOUCHED=TRUE".as_ptr(), ptr::null_mut()] + } else { + [c"ALL_TOUCHED=FALSE".as_ptr(), ptr::null_mut()] + }; + + let geometries_c: Vec<_> = geometries.iter().map(|geo| geo.c_geometry()).collect(); + let burn_values_expanded: Vec<f64> = burn_values + .iter() + .flat_map(|burn| std::iter::repeat_n(burn, bands_i32.len())) + .copied() + .collect(); + + let gt = dataset.geo_transform().map_err(|_e| { + GdalError::BadArgument( + "Missing geotransform: only geotransform-based affine rasterize is supported" + .to_string(), + ) + })?; + let inv_gt = gt.invert().map_err(|_e| { + GdalError::BadArgument( + "Non-invertible geotransform: only geotransform-based affine rasterize is supported" + .to_string(), + ) + })?; + let mut arg = AffineTransformArg { gt, inv_gt }; + + unsafe { + let error = call_gdal_api!( + api, + GDALRasterizeGeometries, + dataset.c_dataset(), + bands_i32.len() as c_int, + bands_i32.as_ptr(), + geometries_c.len() as c_int, + geometries_c.as_ptr(), Review Comment: Added comment to clarify the assumptions on the target platform -- This is an automated message from the Apache Git Service. To respond to the message, please log on to GitHub and use the URL above to go to the specific comment. To unsubscribe, e-mail: [email protected] For queries about this service, please contact Infrastructure at: [email protected]
