RGeo: How to transform a projected coordinate to lat/lon using SRID

Question

I'd like to transform a point in projected coordinates to lat/lon. I've read this blog post about transformations and I came up with the following code:

factory_25831 = RGeo::Cartesian.factory(:srid => 25831)

point_25831 = factory_25831.point(505500.421875, 4699669.4375)

factory_4326 = RGeo::Geographic.spherical_factory(:srid => 4326)

point_4326 = RGeo::Feature.cast(point_25831, :factory => factory_4326, :project => true)

which gives me:

=> #<RGeo::Geographic::SphericalPointImpl:0xe35f6c "POINT (60.421875 90.0)">

However I should have obtained something like long: 3.0668887° / lat: 42.4493345°.

What am I doing wrong?

EDIT:

So apparently I "can" transform a point, but only providing full proj4 and WKT specs, which somehow loses the goodness of providing only the srid code. My code receives a random geometry with a srid as metadata - I therefore cannot hard-code the specification.

proj4_25831 = '+proj=utm +zone=31 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'

wkt_25831 = 'PROJCS["ETRS89 / UTM zone 31N",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","25831"]]'

factory_25831 = RGeo::Cartesian.factory(:srid => 25831,
  :proj4 => proj4_25831, :coord_sys => wkt_25831)



proj4_4326 = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'

wkt_4326 = <<WKT
  GEOGCS["WGS 84",
    DATUM["WGS_1984",
      SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
      AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
      AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]
WKT


factory_4326 = RGeo::Geographic.spherical_factory(:srid => 4326,
  :proj4 => proj4_4326, :coord_sys => wkt_4326)

score 2 · Answer 1 · answered Oct 03 '20 at 15:19

You need to tell the factory where get projection information. This done ay passing a spatial reference database to the factory. With the database, you can look up the dynamically.

srs_database = RGeo::CoordSys::SRSDatabase::Proj4Data.new('epsg', cache: true)

factory_25831 = RGeo::Cartesian.factory(srid: 25831, srs_database: srs_database)

point_25831 = factory_25831.point(505500.421875, 4699669.4375)

factory_4326 = RGeo::Geographic.spherical_factory(srid: 4326, srs_database: srs_database)

point_4326 =  RGeo::Feature.cast(point_25831, factory: factory_4326, project: true)

Coords will be returned accuratly:

point_4326
=> #<RGeo::Geographic::SphericalPointImpl:0x17c "POINT (3.066888712691441 42.44933441408618)">

RGeo: How to transform a projected coordinate to lat/lon using SRID

1 Answers1