proj4rb icon indicating copy to clipboard operation
proj4rb copied to clipboard

Transformation does not work as expected

Open ollietreend opened this issue 10 months ago • 12 comments

I'm trying to convert coordinates from EPSG:27700 to EPSG:4326 using proj4rb.

I have the latest versions of the proj library (installed with brew install proj) and the proj4rb gem installed.

$ proj
Rel. 9.4.0, March 1st, 2024
usage: proj [-bdeEfiIlmorsStTvVwW [args]] [+opt[=arg] ...] [file ...]
$ gem list proj4rb

*** LOCAL GEMS ***

proj4rb (4.1.1)

✅ Using cs2cs command line tool

When I use the command line tool cs2cs provided by proj, I get the expected results:

$ echo "533498 181201" | cs2cs -f '%.10f' EPSG:27700 EPSG:4326
51.5139702631	-0.0775045667 0.0000000000

Easting 533498, Northing 181201 translates to Longitude -0.0775045667, Latitude 51.5139702631. This result is consistent with online converter tools such as the British Geological Survey's Coordinate converter tool.

❌ Using proj4rb

But when I try to perform the same action using proj4rb, I get a completely different result.

require "proj"

transform = Proj::Transformation.new("EPSG:27700", "EPSG:4326")
coordinate = Proj::Coordinate.new(e: 533498, n: 181201)
transformed_coordinate = transform.forward(coordinate)

latitude = transformed_coordinate.y # Latitude is stored in y attribute
longitude = transformed_coordinate.x # Longitude is stored in x attribute

puts "Latitude: #{latitude}"
puts "Longitude: #{longitude}"

When I execute this, I get the following result:

$ ruby proj4.rb
Latitude: 5.0e-324
Longitude: 2.19144499e-314

The returned lat/lng values are effectively zero. If I format them to 10 decimal places (the same as cs2cs output), they are both 0.0000000000.

I need to format them to 400 decimal places to actually see numbers come out:

puts "Latitude: #{"%.400f" % latitude}"
puts "Longitude: #{"%.400f" % longitude}"

which produces:

Latitude: 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000049406564584124654417656879286822137236505980261432476442558568250067550727021
Longitude: 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000219134797539321609555612014012053752743206045798867713541392373366376270284867081861393

What I've tried so far

I wondered if the proj4rb gem was using a different installation of the proj library – perhaps there could have been a broken installation hiding somewhere.

The cs2cs command is located here:

$ realpath $(which cs2cs)
/opt/homebrew/Cellar/proj/9.4.0/bin/cs2cs

And the proj4rb gem seems to be using the same installation path:

require "proj"
Proj::Context.current.database.path
# => "/opt/homebrew/Cellar/proj/9.4.0/share/proj/proj.db"

So from what I can see, they're using the exact same installation. That can't be the issue.

What am I doing wrong?

I'm really confused by this. @cfis do you have any idea what I'm doing wrong?

For what it's worth, I've also noticed that the proj4rb test suite doesn't pass on my machine. I get a whole bunch of failures across pretty much every test file.

ollietreend avatar Apr 04 '24 19:04 ollietreend

Works for me here if I do this:

coordinate = Proj::Coordinate.new(x: 533498, y: 181201)
puts coordinate
v0: 51.513968813644965, v1: -0.0775306317158091, v2: 0.0, v3: 0.0

If I do what you have:

coordinate = Proj::Coordinate.new(e: 533498, n: 181201)
puts coordinate
v0: 49.76680723514262, v1: -7.55715980690519, v2: 0.0, v3: 0.0

The reason is that code requires that you set :e, :n and :u. See https://github.com/cfis/proj4rb/blob/master/lib/proj/coordinate.rb#L62. If not, then no coordinates get set at all. So you are asking to transform an array of values of []. So you get weird numbers.

Really that big if statement should be rewritten where each parameter name (:x, :y, etc) is assigned an index of 0 through 4 and then mapped by index into the Api::PJ_COORD structure. Would be a fairly quick fix if you have any interest in making a PR.

cfis avatar Apr 05 '24 04:04 cfis

@cfis thanks for the quick response! 🙇🏻

I'm afraid I still can't get it to work:

require "proj"

transform = Proj::Transformation.new("EPSG:27700", "EPSG:4326")
coordinate = Proj::Coordinate.new(x: 533498, y: 181201)

puts coordinate
# v0: 533498.0, v1: 181201.0, v2: 0.0, v3: 0.0

puts transform.forward(coordinate)
# v0: 1.7733993367673755e+77, v1: 2.1601183626e-314, v2: 3.0154665634e-314, v3: 1.5e-323

What output do you get when you run the above code?

Based on the code snippets in your previous comment, it looks like you're getting 'good' results whereas I'm still getting something very odd. I'm expecting I should receive meaningful values for v0 and v1, and v2: 0.0, v3: 0.0 because I only provided two params.

ollietreend avatar Apr 05 '24 11:04 ollietreend

In fact, even more confusingly, I consistently get inconsistent results. It seems like the same transformation in a loop always results in the first value being different from the others... 🤔

It feels like something very bizarre is going on, but I'm struggling to know what it could be.

require "proj"

transform = Proj::Transformation.new("EPSG:27700", "EPSG:4326")
coordinate = Proj::Coordinate.new(x: 533498, y: 181201)

5.times { puts transform.forward(coordinate) }
# v0: 5.0e-324, v1: 5.2e-322, v2: 3.025252423e-314, v3: 1.5e-323
# v0: 3.0252521343e-314, v1: 2.1462617777e-314, v2: 3.025252423e-314, v3: 1.5e-323
# v0: 3.0252521343e-314, v1: 2.1462617777e-314, v2: 3.025252423e-314, v3: 1.5e-323
# v0: 3.0252521343e-314, v1: 2.1462617777e-314, v2: 3.025252423e-314, v3: 1.5e-323
# v0: 3.0252521343e-314, v1: 2.1462617777e-314, v2: 3.025252423e-314, v3: 1.5e-323

ollietreend avatar Apr 05 '24 11:04 ollietreend

I've even tried this in Docker to remove the possibility of something weird in my local setup. But I'm still getting the same results 😞

Here is what I've used, and the output I got:

https://gist.github.com/ollietreend/286d2e28806c53c6e022934fbccb0b66

ollietreend avatar Apr 05 '24 12:04 ollietreend

Running the code you provided works for me. I tested on both Windows and Linux (Fedora 39).

Also built and your ran your docker file (thanks for providing!):

docker run proj_test
v0: 533498.0, v1: 181201.0, v2: 0.0, v3: 0.0
v0: 51.513968813644965, v1: -0.0775306317158091, v2: 0.0, v3: 0.0

So sorry, not very helpful...but I have no idea why you are having issues.

Maybe a corrupt proj database? (but then the docker file should avoid that problem).

cfis avatar Apr 05 '24 23:04 cfis

Any luck tracking this down?

cfis avatar Apr 09 '24 07:04 cfis

Sadly, no 😞

The most bizarre thing is that the Docker image doesn't work, and yet it does for you. It makes me wonder if this could be down to the CPU architecture? I'm using a MacBook Pro with M3 processor, so everything's running on ARM. I don't know enough about the internal of this gem to know if that would make a difference or not.

For my particular use case, I need to batch transform about 20k coordinates as part of an ETL process. This runs once a day as a background task, and performance isn't critical.

Since the cs2cs command line tool (provided as part of the PROJ library) works correctly, I've switched to calling that from within my Rails app. I found that I can perform all 20k transformations in a loop using a single cs2cs process by writing to and reading from the IO stream sequentially. This performs surprisingly well for my needs – it does 20k transformations in ~0.2 seconds.

To give you an idea of what I've got working:

class CoordinateTransformer
  # Source coordinates are in British National Grid Easting/Northing format.
  # https://epsg.io/27700
  SOURCE_CRS = "EPSG:27700".freeze

  # We need coordinates in WGS84 format, the standard Latitude/Longitude
  # coordinate system used in GPS and online mapping tools.
  # https://epsg.io/4326
  TARGET_CRS = "EPSG:4326".freeze

  def initialize
    # cs2cs is provided as part of the PROJ library
    # Mac: brew install proj
    # Linux (Debian): apt-get install proj-bin
    # Usage: https://proj.org/en/9.4/apps/cs2cs.html
    @cs2cs = IO.popen(["cs2cs", "-d", "10", SOURCE_CRS, TARGET_CRS], "r+")
  end

  def transform(easting:, northing:)
    @cs2cs.write "#{easting} #{northing}\n"
    output = @cs2cs.gets.split(" ")
    { latitude: output[0], longitude: output[1] }
  end

  def close
    @cs2cs.close
  end
end

As I said, I'm quite impressed by the performance of this approach. I assume it's because it's cheap to write to the existing process rather than continually spawning new ones.

require "benchmark"

Benchmark.measure do
  ct = CoordinateTransformer.new
  20_000.times do |i|
    ct.transform(easting: i, northing: i)
  end
  ct.close
end

# =>   0.048043   0.050908   0.207849 (  0.200632)

It's frustrating I couldn't get the proj4rb gem to work 😢 thank you for your help digging into this with me though @cfis.

ollietreend avatar Apr 09 '24 11:04 ollietreend

Glad you have a solution.

The proj4 gem uses FFI, so does not include any native code (so nothing to compile) and thus is architecture independent. Of course the FFI gem is compiled against the ffi library. But I assume those are well tested and work fine on M3. But I don't have any other explanation for what you seeing.

Unfortunately I don't have a M3 (or M1/M2) computer so can't test for myself.

cfis avatar Apr 09 '24 19:04 cfis

@cfis interestingly, I've just tried running the same Docker image under x86-64/AMD64 and it produces different results compared to ARM 👀

❌ ARM64 (e.g. Apple Silicon) produces bad results

$ docker run --rm --platform linux/arm64 $(docker build -q --platform linux/arm64 .)
v0: 533498.0, v1: 181201.0, v2: 0.0, v3: 0.0
v0: 4.4748e-320, v1: 0.0, v2: 1.0e-323, v3: 1.39067016643389e-309

✅ AMD64 (e.g. Intel) produces good results

docker run --rm --platform linux/amd64 $(docker build -q --platform linux/amd64 .)
v0: 533498.0, v1: 181201.0, v2: 0.0, v3: 0.0
v0: 51.513968813644965, v1: -0.0775306317158091, v2: 0.0, v3: 0.0

To be clear: I'm running both these commands on the same MacBook with M3 processor. The only difference is that I've changed the --platform option.

By default, Docker runs under --platform linux/arm64 on an M3 (ARM) chip. But when you specify --platform linux/amd64 it'll run as AMD64 (presumably using Rosetta).

I'm not entirely sure what this means for proj4rb. I don't know enough about Ruby FFI to know whether this is a bug in proj4rb or FFI itself.

ollietreend avatar Apr 10 '24 13:04 ollietreend

In fact, I think this might just be a general macOS compatibility issue, rather than being specific to the newer Apple Silicon chips. I've forked this repo and added macOS GitHub-hosted runners to the CI workflow, and the test suite fails in the same way it does on my MacBook.

I haven't opened a Pull Request, but you can see what I changed here: https://github.com/cfis/proj4rb/compare/master...ollietreend:proj4rb:macos-runners

You can see the results of this CI workflow here: https://github.com/ollietreend/proj4rb/actions/runs/8634089078

Notably, the macos-13 runner is Intel whereas the macos-14 runner is Apple Silicon (ARM). But both test suites fail with the same errors.

ollietreend avatar Apr 10 '24 15:04 ollietreend

Thanks for keeping at this. I do have an intel Mac, I will give that a try.

That looks like a great update to the test matrix, can you create an MR?

cfis avatar Apr 10 '24 20:04 cfis

@cfis No problem – I've just opened a pull request. See #24.

Although I'm not sure it'll be possible to merge it, since the tests fail on macOS runners 🙈

ollietreend avatar Apr 11 '24 21:04 ollietreend