Reproject from EPSG(4326) to EPSG(29902 seem to give the wrong answer

[AquaPore]

Hi,
A great thanks for this usefull tool.
I believe that there are errors with the reprojection, I am wandering if I am doing anything wrong.

With kindness,
Joseph A.P. POLLACCO

Latitude =[51.629893, 51.62388, 51.621277, 51.622027]

Longitude = [-8.770467, -8.77321, -8.790204, -8.804208]

Path  = raw"D:\TEMP\Gis2\test16_points.shp"

Points = GeoDataFrames.GeoInterface.Point.(Longitude, Latitude, )
Df = GeoDataFrames.DataFrame(customgeomcolumn=Points, name="Datapoints")
Df_Projected = GeometryOps.reproject(Df, EPSG(4326), EPSG(29902))
GeoDataFrames.write(Path, Df_Projected; geometrycolumn=:customgeomcolumn, crs=EPSG(29902))

[asinghvi17]

julia> projected = GO.reproject(Points, EPSG(4326), EPSG(29902))
4-element Vector{Tuple{Float64, Float64}}:
 (146699.21717428684, 42148.195196917804)
 (146502.20607215434, 41481.18262444617)
 (145322.18944846842, 41204.1504572267)
 (144353.2324970336, 41298.170862944826)

julia> GO.reproject(projected, EPSG(29902), EPSG(4326))
4-element Vector{Tuple{Float64, Float64}}:
 (-8.770466994026089, 51.62989299767637)
 (-8.773209994030237, 51.62387999767011)
 (-8.790203994048358, 51.62127699766767)
 (-8.804207994062663, 51.62202699766869)

julia> Points
4-element Vector{Tuple{Float64, Float64}}:
 (-8.770467, 51.629893)
 (-8.77321, 51.62388)
 (-8.790204, 51.621277)
 (-8.804208, 51.622027)

They seem to round trip correctly, at least. What specifically seems to be wrong with them?

[AquaPore]

Thanks, when I save the shapefile the output points are way out when projected in the EPSG:29902 - TM65 / Irish Grid system.

[evetion]

Thanks for making an issue! It seems your custom geometry column is not recognized by reproject. We should add a warning when this happens. You can fix your example by doing:

metadata!(Df, "GEOINTERFACE:geometrycolumns", (:customgeomcolumn,); style = :note))

[asinghvi17]

EDIT: Ah, I didn't see @evetion's comment above before I posted.

Yeah reproject doesn't seem to have a geometrycolumn argument yet - I will fix that in GeometryOps.

Hmm according to https://epsg.io/29902 the axis order is "correct" (east, north) which corresponds to the order we have in Julia (x, y).

If I plot it they also look to be in the same place.

Which viewer are you using that these look bad in?

[AquaPore]

Hi thanks for your help. The following code does not seem to be doing anything:
Df_Projected = GeoDataFrames.reproject(Df, GeoDataFrames.EPSG(29902), GeoDataFrames.EPSG(4326))

I am a bit confused will you kindly let me know how I shoudl write the code:

Many thanks for any help you may provide

using GeoDataFrames

Latitude =[51.6298370, 51.629893, 51.62388, 51.621277,51.622027]
Longitude = [-8.7704106, -8.770467, -8.77321, -8.790204, -8.804208]

Path  = raw"D:\TEMP\Gis2\test17_points.shp"

Points = GeoDataFrames.GeoInterface.Point.(Longitude, Latitude)
Df = GeoDataFrames.DataFrame(customgeomcolumn=Points, name="Datapoints"; crs=EPSG(4326))
GeoDataFrames.metadata!(Df, "GEOINTERFACE:geometrycolumns", (:customgeomcolumn,); style = :note)
Df_Projected = GeoDataFrames.reproject(Df, GeoDataFrames.EPSG(29902), GeoDataFrames.EPSG(4326))
GeoDataFrames.write(Path, Df_Projected; geometrycolumn=:customgeomcolumn, crs=GeoDataFrames.EPSG(29902))

Error code

Can't convert EPSG{1} to an OGRFieldType. Please report an issue.

Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:44
 [2] write(::GeoDataFrames.ArchGDALDriver, fn::String, table::DataFrame; layer_name::String, crs::EPSG{1}, driver::Nothing, options::Dict{String, String}, geom_columns::Nothing, geometrycolumn::Symbol, chunksize::Int64, update::Bool, kwargs::@Kwargs{})
   @ GeoDataFrames C:\Users\jpollacco.local\.julia\packages\GeoDataFrames\TjyE4\src\io.jl:223
 [3] write
   @ C:\Users\jpollacco.local\.julia\packages\GeoDataFrames\TjyE4\src\io.jl:151 [inlined]
 [4] write(driver::GeoDataFrames.ShapefileDriver, fn::String, table::DataFrame; kwargs::@Kwargs{geometrycolumn::Symbol, crs::EPSG{1}})
   @ GeoDataFrames C:\Users\jpollacco.local\.julia\packages\GeoDataFrames\TjyE4\src\io.jl:143
 [5] write
   @ C:\Users\jpollacco.local\.julia\packages\GeoDataFrames\TjyE4\src\io.jl:141 [inlined]
 [6] write(fn::String, table::DataFrame; kwargs::@Kwargs{geometrycolumn::Symbol, crs::EPSG{1}})
   @ GeoDataFrames C:\Users\jpollacco.local\.julia\packages\GeoDataFrames\TjyE4\src\io.jl:133
 [7] top-level scope
   @ d:\JOE\MAIN\MODELS\WFLOW\WflowDataJoe\WflowRaster.jl\src\Other\jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_X12sZmlsZQ==.jl:12

[asinghvi17]

Once you have set the metadata GeometryOps.reproject should work.

[AquaPore]

A great thanks for your help,

I tried to implement your correction, but i was not able to get it working, and I will be ever so gratefull if you could point me to the right direction?

With kindness,
Joseph A.P. Pollacco

using GeoDataFrames

Latitude =[51.6298370, 51.629893, 51.62388, 51.621277,51.622027]
Longitude = [-8.7704106, -8.770467, -8.77321, -8.790204, -8.804208]

Path  = raw"D:\TEMP\Gis2\test18_points.shp"

Points = GeoDataFrames.GeoInterface.Point.(Longitude, Latitude)
Df = GeoDataFrames.DataFrame(customgeomcolumn=Points, name="Datapoints"; crs=EPSG(4326))
Df2 = metadata!(Df, "GEOINTERFACE:geometrycolumns", (:customgeomcolumn,); style = :note)
Df_Projected = GeometryOps.reproject(Df2, GeoDataFrames.EPSG(4326), GeoDataFrames.EPSG(29902))
GeoDataFrames.write(Path, Df_Projected; geometrycolumn=:customgeomcolumn, crs=GeoDataFrames.EPSG(29902)) ```

## Error code
'''julia
Can't onvert EPSG{1} to an OGRFieldType. Please repoty an issue

[evetion]

Your Df has a crs column now due to:
Df = GeoDataFrames.DataFrame(customgeomcolumn=Points, name="Datapoints"; crs=EPSG(4326))

Your snippet simplified:

Df = GeoDataFrames.DataFrame(customgeomcolumn=Points, name="Datapoints")
Df2 = metadata!(Df, "GEOINTERFACE:geometrycolumns", (:customgeomcolumn,); style = :note) # required because of the custom geometry column name
Df_Projected = GeometryOps.reproject(Df2, GeoDataFrames.EPSG(4326), GeoDataFrames.EPSG(29902)) # this set the crs metadata
GeoDataFrames.write(Path, Df_Projected)

[AquaPore]

A great thanks for your kindness and help @evetion and @asinghvi17,

It works now beautifully.

I just noticed that your convention is Point.(Latitude_Y, Longitude_X) and not Point.(Latitude_X, Longitude_Y)

Shall I close the case?

With kindness,
Joseph

using GeoDataFrames

Latitude_Y =[51.629893
51.62388
51.62388
51.621277
51.622027
]
Longitude_X = [-8.770467
-8.77321
-8.77321
-8.790204
-8.804208
]
Site =["5_ASU1"
"5_ASU2"
"5_ASU2"
"5_ASU3"
"5_ASU4"
]
Path  = raw"D:\TEMP\Gis2\test22_points.shp"

Points = GeoDataFrames.GeoInterface.Point.(Latitude_Y, Longitude_X)

Df = GeoDataFrames.DataFrame(Coordinates=Points, Site=Site)
Df = GeoDataFrames.metadata!(Df, "GEOINTERFACE:geometrycolumns", (:Coordinates,); style = :note) # required because of the custom geometry column name
Df_Projected = GeometryOps.reproject(Df, GeoDataFrames.EPSG(4326), GeoDataFrames.EPSG(29902)) # this set the crs metadata
GeoDataFrames.write(Path, Df_Projected)

[evetion]

I just noticed that your convention is Point.(Latitude_Y, Longitude_X) and not Point.(Latitude_X, Longitude_Y)

I'm not sure what you mean here with latitude_y/latitude_x, but I think your question is how to order the coordinates. That depends on the CRS(!), per https://proj.org/en/stable/faq.html#why-is-the-axis-ordering-in-proj-not-consistent. For geodetic stuff, it's often latitude, longitude, whereas most projected things use x, y. It can be confusing.

edit: Note that there's a always_xy (false by default) keyword that you can pass to force the ordering (i.e. longitude, latitude) for your input.

[asinghvi17]

We set always_xy to true in GeometryOps by default. So the ordering is always

(long, lat)

or

(x, y)

[AquaPore]

Dear @evetion and @asinghvi17 ,

A great thanks for helping to remove the wild guess of the convention used which indeed is confusing.

I am still a bit confused. Therefore if the convention is as pointed out by asinghvi (long, lat) and (x, y) then I guess that an option is needed in the Point. Currently with the EPSG_Input=4326 has the following convention:

GeoDataFrames.GeoInterface.Point.(Latitude_Y, Longitude_X)

May be there should be an option in Point, always_xy=true such that

GeoDataFrames.GeoInterface.Point.( Longitude_X, Latitude_Y, always_xy=true)

This is just a taught,

With kindness,
J.A.P. Pollacco

[evetion]

Hi @AquaPore, thanks for the ideas. I just merged a PR here that will set always_xy=true by default for GeoDataFrames as well. It seems good to unify long/lat | x/y accross the JuliaGeo ecosystem, so we have no confusion about it. That means you should not do Point(latitude, longitude).

FWIW, you can use kwargs:

julia> latitude, longitude = 52, 4
(52, 4)

julia> GeoDataFrames.GeoInterface.Point(Y=latitude, X=longitude)
Point{false, false}((4, 52))