rgeo / rgeo-proj4

Proj4 extension for rgeo.
MIT License
13 stars 14 forks source link

Regression in projection definition #16

Closed inkstak closed 1 year ago

inkstak commented 3 years ago

There is a regression between rgeo-proj4 2.0.1 and 3.0.1

The following example works in 2.0.1:

#  rgeo-proj4 2.0.1

source_projection = RGeo::CoordSys::Proj4.new('+proj=longlat +datum=WGS84 +no_defs')
target_projection = RGeo::CoordSys::Proj4.new('+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ')

RGeo::CoordSys::Proj4.transform_coords(source_projection, target_projection, 2.3347587142445168, 48.738770796515, nil)
# => [651078.9599999869, 6848945.02999845]

But not in 3.0.1:

#  rgeo-proj4 3.0.1

source_projection = RGeo::CoordSys::Proj4.new('+proj=longlat +datum=WGS84 +no_defs')
target_projection = RGeo::CoordSys::Proj4.new('+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs')

RGeo::CoordSys::Proj4.transform_coords(source_projection, target_projection, 2.3347587142445168, 48.738770796515, nil)
# proj_create_operations: source_crs is not a CRS
# => nil

To make it work, I had to add +type=crs to both projection or use projection name (which is better)

#  rgeo-proj4 3.0.1

source_projection = RGeo::CoordSys::Proj4.new('+proj=longlat +datum=WGS84 +no_defs +type=crs')
target_projection = RGeo::CoordSys::Proj4.new('+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs')

RGeo::CoordSys::Proj4.transform_coords(source_projection, target_projection, 2.3347587142445168, 48.738770796515, nil)
# => [651078.9599999869, 6848945.02999845]

source_projection = RGeo::CoordSys::Proj4.new('EPSG:4326')
target_projection = RGeo::CoordSys::Proj4.new('EPSG:2154')

RGeo::CoordSys::Proj4.transform_coords(source_projection, target_projection, 2.3347587142445168, 48.738770796515, nil)
# => [651078.9599999869, 6848945.02999845]

FYI, epsg.io doesn't provide definition with +type=crs: https://epsg.io/2154.proj4 (spatialreference.org is down, so I cannot double check)

It may be a issue with proj4 and doing so should be noted as a breaking change ?

inkstak commented 3 years ago

BTW, RGeo still uses few projections that cannot be projected with this new version:

https://github.com/rgeo/rgeo/blob/master/lib/rgeo/geographic/interface.rb#L464

source_projection = RGeo::Geographic.simple_mercator_factory.proj4
# => #<RGeo::CoordSys::Proj4:0x16b34 "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs">
keithdoggett commented 3 years ago

Hi, @inkstak thanks for the detailed issue.

In regards to the first part about +type=crs, yes that's necessary and I can add a note to the README. PROJ has been updated a lot since v4, so the definitions have changed and in newer versions, you need to add +type=crs to define a CRS. I was considering automatically adding that on creation, but thought it might be confusing since it's hiding changes in the definition from the user. I could potentially revisit that, but I'm still don't think it's the best idea and could constrain us down the line if we wanted to add non-CRS functionality.

FWIW PROJ no longer recommends the use of proj strings, so ideally I'd migrate away from them entirely. Here are their notes on CRS creation: https://proj.org/development/reference/functions.html#c.proj_create.

Thanks for pointing out the proj strings used in the core library. I can update those definitions soon.

Also, the simple_mercator_factory should still work even though it's proj4 attribute is wrong since the projection is implemented in Ruby. But I'll update it's definition to be correct.

inkstak commented 3 years ago

Ok, I better understand the issue.

We are already replacing proj strings by their EPSG equivalents whenever it's possible, but this is not always the case when it comes to deal with user files.

For example : we use GDAL OGR drivers - https://gdal.org/drivers/vector/index.html - to import several data formats (mostly Shapefile, EDIGEO, ...) via the gdal gem - https://github.com/zhm/gdal-ruby - and the projection is always returned as a proj string. So I hope you'll maintain them as long as possible.

Our workaround is, as you suggest, to check and add +type=crs at the end of the string when it's missing.

keithdoggett commented 3 years ago

Yeah, I think the best thing, for now, is to do that check on your end and add them. I'll update the core gem to use the modern syntax and add a note to the README here.

Just to clarify, as long as PROJ supports proj strings this gem will support them, I meant that I'll replace the hardcoded proj strings with the EPSG codes throughout the core gem.

keithdoggett commented 1 year ago

This has been fixed with the RGeo 3 release and RGeo-Proj4 V4 release