Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add transform_bounds to wrap proj_trans_bounds #155

Merged
merged 3 commits into from
Jun 9, 2023
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
88 changes: 84 additions & 4 deletions src/proj.rs
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,8 @@ use proj_sys::{
proj_context_is_network_enabled, proj_context_set_search_paths, proj_context_set_url_endpoint,
proj_create, proj_create_crs_to_crs, proj_destroy, proj_errno_string, proj_get_area_of_use,
proj_grid_cache_set_enable, proj_info, proj_normalize_for_visualization, proj_pj_info,
proj_trans, proj_trans_array, PJconsts, PJ_AREA, PJ_CONTEXT, PJ_COORD, PJ_DIRECTION_PJ_FWD,
PJ_DIRECTION_PJ_INV, PJ_INFO, PJ_LPZT, PJ_XYZT,
proj_trans, proj_trans_array, proj_trans_bounds, PJconsts, PJ_AREA, PJ_CONTEXT, PJ_COORD,
PJ_DIRECTION_PJ_FWD, PJ_DIRECTION_PJ_INV, PJ_INFO, PJ_LPZT, PJ_XYZT,
};
use std::{
convert, ffi,
Expand Down Expand Up @@ -426,7 +426,6 @@ impl ProjBuilder {
///
///```rust
/// # use approx::assert_relative_eq;
/// extern crate proj;
/// use proj::{Proj, Coord};
///
/// let from = "EPSG:2230";
Expand Down Expand Up @@ -787,7 +786,6 @@ impl Proj {
///
/// ```rust
/// # use approx::assert_relative_eq;
/// extern crate proj;
/// use proj::{Proj, Coord};
///
/// let from = "EPSG:2230";
Expand Down Expand Up @@ -916,6 +914,88 @@ impl Proj {
self.array_general(points, Transformation::Projection, inverse)
}

/// Transform boundary densifying the edges to account for nonlinear transformations along
/// these edges and extracting the outermost bounds.
///
/// Input and output CRS may be specified in two ways:
/// 1. Using the PROJ `pipeline` operator. This method makes use of the [`pipeline`](http://proj4.org/operations/pipeline.html)
/// functionality available since `PROJ` 5.
/// This has the advantage of being able to chain an arbitrary combination of projection, conversion,
/// and transformation steps, allowing for extremely complex operations ([`new`](#method.new))
/// 2. Using EPSG codes or `PROJ` strings to define input and output CRS ([`new_known_crs`](#method.new_known_crs))
///
/// The `densify_pts` parameter describes the number of points to add to each edge to account
/// for nonlinear edges produced by the transform process. Large numbers will produce worse
/// performance.
///
/// The following example converts from NAD83 US Survey Feet (EPSG 2230) to NAD83 Metres (EPSG 26946)
///
/// ```rust
/// # use approx::assert_relative_eq;
/// use proj::{Proj, Coord};
///
/// let from = "EPSG:2230";
/// let to = "EPSG:26946";
/// let ft_to_m = Proj::new_known_crs(&from, &to, None).unwrap();
/// let result = ft_to_m
/// .transform_bounds(4760096.421921, 3744293.729449, 4760196.421921, 3744393.729449, 21)
/// .unwrap();
/// assert_relative_eq!(result[0] as f64, 1450880.29, epsilon=1e-2);
/// assert_relative_eq!(result[1] as f64, 1141263.01, epsilon=1e-2);
/// assert_relative_eq!(result[2] as f64, 1450910.77, epsilon=1e-2);
/// assert_relative_eq!(result[3] as f64, 1141293.49, epsilon=1e-2);
/// ```
///
/// # Safety
/// This method contains unsafe code.
pub fn transform_bounds<F>(
&self,
left: F,
bottom: F,
right: F,
top: F,
densify_pts: i32,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A doc about densify pts would be nice.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added a sentence taken from the pyproj docs

) -> Result<[F; 4], ProjError>
Copy link
Member

@michaelkirk michaelkirk Mar 22, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Though it's pretty much just like pyproj does, it seems a little strange to me that the input and output formats are different. (left, bottom, right, top) -> [F; 4]

We could do like we did with transform and accept some kind of Bounds trait.

Or we could output a bespoke Bounds { left, bottom, right, top }struct (or maybe reuse/consolidate with the existing proj::Area).

One thing I considered, if we did go the Bounds trait route, it'd be natural to want to do:

impl proj::Bounds for geo_types::Rect {
...
}

... but I think that might give us some incorrect results when we span the prime meridian - with how you have it now, it's no problem if left > right. I'm not actually sure what happens if you try it though.

Or we can just leave it. Apparently it's good enough for pyproj.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it seems a little strange to me that the input and output formats are different

That's fair. Open to change it to return a tuple of F or use a struct. I'm not too familiar with existing crate conventions.

when we span the prime meridian

It looks like proj allows the max value to be less than the min when it spans the antimeridian.

https://github.com/OSGeo/PROJ/blob/04615b8dbe9c6da65e6898255adeda16dd907ddc/src/4D_api.cpp#L1536-L1537

where
F: CoordinateType,
{
let mut new_left = f64::default();
let mut new_bottom = f64::default();
let mut new_right = f64::default();
let mut new_top = f64::default();
let err;

unsafe {
proj_errno_reset(self.c_proj);
let _success = proj_trans_bounds(
self.ctx,
self.c_proj,
PJ_DIRECTION_PJ_FWD,
left.to_f64().ok_or(ProjError::FloatConversion)?,
bottom.to_f64().ok_or(ProjError::FloatConversion)?,
right.to_f64().ok_or(ProjError::FloatConversion)?,
top.to_f64().ok_or(ProjError::FloatConversion)?,
&mut new_left,
&mut new_bottom,
&mut new_right,
&mut new_top,
densify_pts,
);
err = proj_errno(self.c_proj);
}

if err == 0 {
Ok([
F::from(new_left).ok_or(ProjError::FloatConversion)?,
F::from(new_bottom).ok_or(ProjError::FloatConversion)?,
F::from(new_right).ok_or(ProjError::FloatConversion)?,
F::from(new_top).ok_or(ProjError::FloatConversion)?,
])
} else {
Err(ProjError::Conversion(error_message(err)?))
}
}

// array conversion and projection logic is almost identical;
// transform points in input array into PJ_COORD, transform them, error-check, then re-fill
// input slice with points. Only the actual transformation ops vary slightly.
Expand Down