DHI / mikecore-python

MIKE Core in Python
5 stars 0 forks source link

Projected to Geographical conversion breaks with Lambert_Conformal_Conic_2SP #26

Open stephankistner opened 1 year ago

stephankistner commented 1 year ago

Hi,

I had issues in mikeio when trying to create a Dfs2 file with a _Lambert_Conformal_Conic2SP projection. I can do this succesfully when trying manually (i.e. in MIKE Zero GUI), and in mikecore. So I traced the issue to the mikecore.projections.MapProjection class used to setup the Dfs2 projection info for the Dfs2 builder: when it tries to convert projected coordinates it runs into numerical issues. See example below:

Fail case: Lambert_Conformal_Conic_2SP projection

from mikecore.Projections import MapProjection

# The projected and geographic coordinates here are the same
# This can be converted back and forth using the Datum Converter tool in MIKE Zero
origin_projected = [  -65900.32888812723, -678217.039068304]
origin_geographic = [  24.01529,-34.987278]

# Custom projection string
proj_str = 'PROJCS["Custom",GEOGCS["Unused",DATUM["User defined",SPHEROID["WGS 1984",6378137,298.257223563]],'+\
    'PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],'+\
    'PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",24.751587],'+\
    'PARAMETER["Standard_Parallel_1",-40],PARAMETER["Standard_Parallel_2",-10],'+\
    'PARAMETER["Latitude_Of_Origin",-28.702538],UNIT["Meter",1]]'

# Try conversion using MIKE MapProjection
proj_core = MapProjection.Create(proj_str)

# Projected to Geographic Fails
print(proj_core.Proj2Geo(*origin_projected))

# Geographic to Projected works
print(proj_core.Geo2Proj(*origin_geographic))
>>> (444.87039857071056, nan)
>>> (-65900.2668223073, -678217.0403999202)

As you can see the _Lambert_Conformal_Conic2SP coordinate pairs fail when trying to converted from projected to geographical coordinates, while the reverse works just fine. I haven't tested this in the .NET/C# code, but I'm confident this isn't correct behaviour.

Working case: Transverse_Mercator projection

from mikecore.Projections import MapProjection

# Coordinate info
# The projected and geographic coordinates here are the same
# This can be converted back and forth using the Datum Converter tool in MIKE Zero
origin_projected = [  -89907.8133323,-3873624.5460040]
origin_geographic = [  24.01529,-34.987278]

# Custom projection string
proj_str = 'PROJCS["WG25",GEOGCS["Unused",DATUM["User defined",SPHEROID["WGS 1984",6378137,298.257223563]],'+\
    'PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],'+\
    'PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",25],PARAMETER["Scale_Factor",1],'+\
    'PARAMETER["Latitude_Of_Origin",0],UNIT["Meter",1]]'

# Try conversion using MIKE MapProjection
proj_core = MapProjection.Create(proj_str)

# Projected to Geographic Fails
print(proj_core.Proj2Geo(*origin_projected))

# Geographic to Projected works
print(proj_core.Geo2Proj(*origin_geographic))
>>> (24.015290000000082, -34.9872779999965)
>>> (-89907.81333230426, -3873624.5460043526)

Success - works as expected!

stephankistner commented 1 year ago

Solution/Fix:
The problem is the MzCart.dll version 20.0.0.15310. It's encrypted, so I can't trace the issue further, but I replaced it with the version 21.0.0.16308 from the MIKE 2023 SDK and that solves the issue.

JesperGr commented 1 year ago

Thanks for finding a solution. I think it is time to bump the references to the MIKE Core/SDK components from 20.0 to 21.0. I will put that on my todo-list.