2

I'm trying to use those parameters and I receive a difference of ~0.5m

ALON and BSHM APN Points:

coords = [
  (dms2rad([31,42,28.5694]),  dms2rad([34,36,23.92196]),  55.671),
  (dms2rad([32,46,44.34472]), dms2rad([35,1,22.74061]), 225.046),
]

Transform:

def lla_to_ig05(lon, lat, alt):
    import pyproj
    lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
    ig05 = pyproj.Proj("+proj=tmerc +lat_0=31.7343936111111 +lon_0=35.2045169444445 +k=1.0000067 +x_0=219529.584 +y_0=626907.39 +ellps=GRS80 +towgs84=-24.002400,-17.103200,-17.844400,-0.33077,-1.852690,1.669690,5.424800 +units=m +no_defs")
    x, y, h = pyproj.transform(lla, ig05, lon, lat, alt, radians=True)
    return x, y, h
for lat, lon, alt in coords:
  print(pd.DataFrame(lla_to_ig05(lon, lat, alt))) # pd for pandas
  print("")

Output:

               0
0  162790.443392
1  624087.953905
2      55.671000

               0
0  202456.903152
1  742717.095893
2     225.046000

Real numbers:

               0
0  162789.564
1  624088.355
2      55.671000

               0
0  202456.453
1  742718.779
2     225.046000

Any hint will be very helpful!

2
  • Hi Sasha, welcome to GIS StackExchange. Please post the question here instead of just linking to an answer with a question in it. You really should delete the answer, as it is not an answer. You can still link to a question if it is related. But your question should contain all the information needed to help people help you. Oct 31, 2019 at 14:25
  • Try changing the signs on the rotations--the 4-6 parameters of the transformation. EPSG lists these with coordinate frame method, but PROJ uses position vector.
    – mkennedy
    Oct 31, 2019 at 20:57

1 Answer 1

0

Using the official EPSG codes provides the correct results (pyproj version: 2.4.0 [PROJ version: 6.2.0]):


>>> from pyproj import CRS. Transformer
>>> cc = CRS("IGD05/12")
>>> cc
<Geographic 2D CRS: EPSG:7139>
Name: IGD05/12
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: Israel
- bounds: (32.99, 29.45, 35.69, 33.53)
Datum: Israeli Geodetic Datum 2005(2012)
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

>>> cc2 = CRS("EPSG:6991")
>>> cc2
<Projected CRS: EPSG:6991>
Name: Israeli Grid 05/12
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Asia - Middle East - Israel and Palestine Territory onshore
- bounds: (34.17, 29.45, 35.69, 33.28)
Coordinate Operation:
- name: Israeli TM
- method: Transverse Mercator
Datum: IG05/12 Intermediate Datum
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

>>> trans = Transformer.from_crs(cc, cc2, always_xy=True)
>>> trans
<Concatenated Operation Transformer: pipeline>
Description: axis order change (2D) + IGD05/12 to  IG05/12 Intermediate CRS + Israeli TM
Area of Use:
- name: Asia - Middle East - Israel and Palestine Territory onshore
- bounds: (34.17, 29.45, 35.69, 33.28)
>>> trans.transform(lon, lat, alt)
(162789.56376679792, 624088.3548852643, 55.671)

2
  • Thanks! that way worked.
    – Sasha
    Nov 7, 2019 at 13:24
  • Glad to hear that it worked 👍 - just got a notification pointing to this question again today and realized that I didn't respond.
    – snowman2
    Apr 12, 2020 at 19:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.