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. – HeikkiVesanto Oct 31 '19 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 '19 at 20:57
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 '19 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 '20 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.