WGS84 <--> UTM projection header-only API #1191

harrism · 2023-06-08T11:17:17Z

Depends on #1132
Closes #1139

Description

Implements WGS to UTM projection (either direction) in CUDA C++. This is based on the code from PROJ, and that library is used for comparison in the tests.

To reviewers, I would start by looking over the test: wgs_to_utm_test.cu. This uses projection_factories.hpp, so review that next, then look at projection.cuh and projection_parameters.cuh. The projection pipeline is applied using a thrust::transform() where the operator is a pipeline made up of a sequence of operations (operation.cuh). A pipeline can be applied forward or inverse. The operations for UTM projection are all in include/cuproj/operation.

Note that the dispatch of operations is a little clunky (switch statement in a device function). This is due to limitations of creating arrays of objects on the host for virtual dispatch on the device. This bears more experimentation in the future to come up with a better approach, but this works fine for now.

~~TODO:~~

Design a way to build pipelines of device operations dynamically and run on GPU
Documentation
More testing

Benchmark and Python API are left to follow-up PRs.

Checklist

I am familiar with the Contributing Guidelines.
New or existing tests cover these changes.
The documentation is up to date with these changes.

jarmak-nv

Very much a novice in the C++/CUDA world - a few comments, but it looks very nice! From the GIS perspective and someone who has used PyProj often, it looks like the backend I'd expect to see 🌐

I'll file an issue to expand the readme once we figure out exactly what will be making it into 23.08 (column API, Python API etc.)

This might be me failing to follow along how this works, so maybe a bunch of nonsense below:

Looking at wgs_to_utm_test.cu we talk about (lat, lon) pairs as the coordinate to be transformed, but as far as I can see coordinates are (x, y) formatted where (lat, lon) is (y, x) formatted (never enjoyed this little "feature" of lat/lon). I think we just need to update any mention of (lat, lon) to (lon, lat).

jarmak-nv · 2023-06-29T01:53:19Z

cpp/cuproj/tests/wgs_to_utm_test.cu

+  // generate (lat, lon) points on a grid between -60 and 60 degrees longitude and
+  // -40 and 80 degrees latitude
+  int num_points_x         = 100;
+  int num_points_y         = 100;
+  coordinate<T> min_corner = {-26.5, -152.5};
+  coordinate<T> max_corner = {-25.5, -153.5};


I think the comment and the generated coord ranges disagree here

Fixed the comment. Regarding lon, lat vs. lat, lon. This is a good point. Currently I assume the data is lat, lon which is the opposite of assumptions in cuSpatial. I inject an axis_swap to convert to lon, lat, in fact, which is what PROJ does. But when the user has data that is lon, lat, that axis_swap is not correct, so this needs to be made a configuration option.

harrism · 2023-06-29T02:13:11Z

I'll file an issue to expand the readme once we figure out exactly what will be making it into 23.08 (column API, Python API etc.)

I was thinking about this -- do we want a column API? Or do we just want to deal with generic Python arrays of some sort? Decision on this will affect what the layer that provides Python binding points will look like.

jarmak-nv · 2023-06-29T02:27:34Z

I was thinking about this -- do we want a column API? Or do we just want to deal with generic Python arrays of some sort? Decision on this will affect what the layer that provides Python binding points will look like.

A good team question! When considering from the cuSpatial point-of-view I think a column API seems right. When considering cuProj alone - that's more interesting. A few questions, ordered in what I think is importance, since fairly unfamiliar:

Which one is more straightforward to build?
Does either approach lend itself better to this solution the cuCIM team posted to accelerate geotiffs using cuPy??
Does using the generic Python arrays pose more work on the cuSpatial side to connect?
Am I right in assuming the generic Python approach would be easier to maybe someday fit with Proj/PyProj?

isVoid · 2023-06-30T19:06:23Z

@harrism when you say generic python array - do you mean arrays that conforms to __cuda_array_interface__? Or it's not limited to that? Just want to clarify since the implementation could vary by a lot.

isVoid

Awesome. This is looking very promising. I'm not pressing request changes because I'm out next week. Don't want to be the blocker of merging.

cpp/cuproj/include/cuproj/operation/clamp_angular_coordinates.cuh

cpp/cuproj/include/cuproj/operation/transverse_mercator.cuh

cpp/cuproj/include/cuproj/projection_factories.hpp

cpp/cuproj/tests/wgs_to_utm_test.cu

isVoid · 2023-07-12T08:53:57Z

cpp/cuproj/include/cuproj/operation/degrees_to_radians.cuh

@@ -57,7 +55,7 @@ class degrees_to_radians : operation<Coordinate> {
  __host__ __device__ Coordinate forward(Coordinate const& coord) const
  {
    using T = typename Coordinate::value_type;
-    return Coordinate{static_cast<T>(coord.x * DEG_TO_RAD), static_cast<T>(coord.y * DEG_TO_RAD)};
+    return Coordinate{coord.x * DEG_TO_RAD<T>, coord.y * DEG_TO_RAD<T>};


I assume we expect other coordinate type than just cuspatial::vec_2d? Otherwise I thought we can just do coord * DEG_TO_RAD<T>.

Yeah, I was concerned about that not working for POD data types, for example. But perhaps I shouldn't worry and just assume proper coordinate types. We can't even use it directly with PJ_COORD because it's a union... https://proj.org/en/9.2/development/reference/datatypes.html#c.PJ_COORD

We could probably provide an adaptor from PJ_COORD in the future.

harrism · 2023-07-20T01:08:58Z

/merge

Closes #1215 Depends on #1191 New benchmark for C++ API for WGS84->UTM transform. Below results were run on a single CPU core for Proj (ignore gbench's claims to the contrary), and on a single H100 80GB GPU. Machine was a DGX H100. ``` (rapids) coder ➜ ~/cuspatial $ CUDA_VISIBLE_DEVICES=1 cpp/build/latest/cuproj/benchmarks/WGS_TO_UTM_BENCH 2023-08-02T00:08:52+00:00 Running cpp/build/latest/cuproj/benchmarks/WGS_TO_UTM_BENCH Run on (224 X 3800 MHz CPU s) CPU Caches: L1 Data 48 KiB (x112) L1 Instruction 32 KiB (x112) L2 Unified 2048 KiB (x112) L3 Unified 107520 KiB (x2) Load Average: 1.70, 6.21, 14.43 --------------------------------------------------------------------------------------------------------------------- Benchmark Time CPU Iterations UserCounters... --------------------------------------------------------------------------------------------------------------------- proj_utm_benchmark/forward_double/100 0.013 ms 0.013 ms 57756 items_per_second=7.95314M/s proj_utm_benchmark/forward_double/1000 0.114 ms 0.114 ms 6118 items_per_second=8.73369M/s proj_utm_benchmark/forward_double/10000 1.24 ms 1.24 ms 588 items_per_second=8.07697M/s proj_utm_benchmark/forward_double/100000 12.0 ms 12.0 ms 58 items_per_second=8.35593M/s proj_utm_benchmark/forward_double/1000000 120 ms 120 ms 6 items_per_second=8.36301M/s proj_utm_benchmark/forward_double/10000000 1213 ms 1213 ms 1 items_per_second=8.24563M/s proj_utm_benchmark/forward_double/100000000 11977 ms 11976 ms 1 items_per_second=8.35038M/s proj_utm_benchmark/forward_double/1000000000 119680 ms 119677 ms 1 items_per_second=8.35579M/s cuproj_utm_benchmark/forward_float/100/manual_time 0.011 ms 0.040 ms 66283 items_per_second=9.47605M/s cuproj_utm_benchmark/forward_float/1000/manual_time 0.012 ms 0.041 ms 56799 items_per_second=81.2019M/s cuproj_utm_benchmark/forward_float/10000/manual_time 0.013 ms 0.042 ms 55571 items_per_second=793.482M/s cuproj_utm_benchmark/forward_float/100000/manual_time 0.013 ms 0.042 ms 53048 items_per_second=7.5779G/s cuproj_utm_benchmark/forward_float/1000000/manual_time 0.027 ms 0.056 ms 25842 items_per_second=36.9063G/s cuproj_utm_benchmark/forward_float/10000000/manual_time 0.170 ms 0.198 ms 4130 items_per_second=58.8554G/s cuproj_utm_benchmark/forward_float/100000000/manual_time 1.60 ms 1.62 ms 439 items_per_second=62.6581G/s cuproj_utm_benchmark/forward_float/1000000000/manual_time 15.9 ms 15.9 ms 44 items_per_second=63.0518G/s cuproj_utm_benchmark/forward_double/100/manual_time 0.012 ms 0.041 ms 57960 items_per_second=8.30297M/s cuproj_utm_benchmark/forward_double/1000/manual_time 0.015 ms 0.044 ms 47605 items_per_second=68.0791M/s cuproj_utm_benchmark/forward_double/10000/manual_time 0.015 ms 0.044 ms 47353 items_per_second=676.864M/s cuproj_utm_benchmark/forward_double/100000/manual_time 0.016 ms 0.045 ms 43394 items_per_second=6.19684G/s cuproj_utm_benchmark/forward_double/1000000/manual_time 0.042 ms 0.070 ms 16621 items_per_second=23.7599G/s cuproj_utm_benchmark/forward_double/10000000/manual_time 0.304 ms 0.332 ms 2302 items_per_second=32.863G/s cuproj_utm_benchmark/forward_double/100000000/manual_time 2.93 ms 2.96 ms 240 items_per_second=34.087G/s cuproj_utm_benchmark/forward_double/1000000000/manual_time 29.3 ms 29.3 ms 24 items_per_second=34.13G/s Authors: - Mark Harris (https://github.com/harrism) Approvers: - Michael Wang (https://github.com/isVoid) URL: #1216

harrism and others added 30 commits May 9, 2023 04:56

Add cmake infrastructure and placeholder files for projection library

cfd69ad

Merge branch 'branch-23.06' into feature/projection/cmake

dcb23ef

Rename to cuproj

7b2e9b9

Merge branch 'branch-23.06' into feature/projection/cmake

f469cf0

Delete cugeo dir

68acf29

show more digits in vec_2d_matcher

c16e735

Add copy of PROJ license

818a3e5

Initial tmerc implementation and test

c1550cd

Merge branch 'branch-23.08' into feature/projection/cmake

c2784a9

Remove _ctk_static_suffix

d8d851d

if TARGET conda_env update

dade592

Merge branch 'branch-23.08' into feature/projection/cmake

10b0590

Update version to 23.08 and PUBLIC -> PRIVATE

88a3760

Merge branch 'feature/projection/cmake' into feature/projection/wgs-utm

c36c7eb

Move PROJ license into the file that uses its code.

2be7746

Remove placeholder CPU code

87172f0

Rename test.

8cdc8dc

Move cuproj into cpp and build it when cuspatial builds.

48ead2e

re-add cuml

b889452

Add sqlite dependency

2073b1f

Merge branch 'branch-23.08' into feature/projection/cmake

d366ad5

Add proj and sqlite dependencies to conda meta.yaml

d59cdc4

Merge branch 'branch-23.08' into feature/projection/cmake

3eba34e

Fix comment

cd68e4e

rename to get_cudf.cmake

5243a5d

Remove GNUInstallDirs

e27f7bb

Copyright

a0ad375

cuspatial-exports --> cuproj-exports

e2a28ea

PROJ namespace

5d395c4

Remove FIND_PACKAGE_ARGUMENTS

517e65b

harrism added feature request New feature or request and removed libcuspatial Relates to the cuSpatial C++ library labels Jun 28, 2023

harrism changed the title ~~WGS84 -> UTM projection header-only API~~ WGS84 <--> UTM projection header-only API Jun 28, 2023

harrism requested review from thomcom and jarmak-nv June 28, 2023 13:38

jarmak-nv reviewed Jun 29, 2023

View reviewed changes

harrism mentioned this pull request Jun 29, 2023

[DOC]: cuProj Readme #1206

Closed

Improve reverse construction logic and fix failing tests

19e9699

github-actions bot added the libcuspatial Relates to the cuSpatial C++ library label Jun 29, 2023

isVoid reviewed Jul 1, 2023

View reviewed changes

harrism added 9 commits July 4, 2023 03:38

Factor out wrapping to -PI..PI and constants

292f85e

Remove commented out code

c23db51

Factor out central meridian/parallel, remove unused comments

ce69abc

Use HOST_DEVICE_EXPECTS and comment instead of #if 0

38fe73c

Improve mismatch printing

cc5afb2

Test multiple places around the globe.

e9471a0

Merge branch 'branch-23.08' into feature/projection/wgs-utm

6e02100

Simplify test vector conversion

2fe01f7

Move grid_generator to cuproj_test

09d0b63

harrism requested a review from isVoid July 11, 2023 20:28

harrism mentioned this pull request Jul 12, 2023

WGS84 to UTM benchmark (C++) #1216

Merged

isVoid approved these changes Jul 12, 2023

View reviewed changes

thomcom approved these changes Jul 19, 2023

View reviewed changes

Merge branch 'branch-23.08' into feature/projection/wgs-utm

feba258

rapids-bot bot merged commit b0657dd into rapidsai:branch-23.08 Jul 20, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

WGS84 <--> UTM projection header-only API #1191

WGS84 <--> UTM projection header-only API #1191

harrism commented Jun 8, 2023 •

edited

Loading

jarmak-nv left a comment

jarmak-nv Jun 29, 2023

harrism Jun 29, 2023

harrism commented Jun 29, 2023

jarmak-nv commented Jun 29, 2023

isVoid commented Jun 30, 2023

isVoid left a comment

isVoid Jul 12, 2023

harrism Jul 12, 2023

harrism Jul 20, 2023

harrism commented Jul 20, 2023

WGS84 <--> UTM projection header-only API #1191

WGS84 <--> UTM projection header-only API #1191

Conversation

harrism commented Jun 8, 2023 • edited Loading

Description

Checklist

jarmak-nv left a comment

Choose a reason for hiding this comment

jarmak-nv Jun 29, 2023

Choose a reason for hiding this comment

harrism Jun 29, 2023

Choose a reason for hiding this comment

harrism commented Jun 29, 2023

jarmak-nv commented Jun 29, 2023

isVoid commented Jun 30, 2023

isVoid left a comment

Choose a reason for hiding this comment

isVoid Jul 12, 2023

Choose a reason for hiding this comment

harrism Jul 12, 2023

Choose a reason for hiding this comment

harrism Jul 20, 2023

Choose a reason for hiding this comment

harrism commented Jul 20, 2023

harrism commented Jun 8, 2023 •

edited

Loading