Stereographic projection of WGS84 ellipsoid on a plane[python]

Question

I would like to ask you for help in following task.
What I need to do is to derive local x,y coordinates from latitude, longitude, applying stereographic projection of WGS84 ellipsoid on a plane, when the following parameters are known:

scale = 0.999790760649094
x0 = 41614.2107651061
y0 = 17150.1692507701
L0 = 18.4514121640352
B0 = 54.4353877827032

There's also a Q parameter but at the moment I hava no idea what is and what for it is:

Q = -0.0203604068846608

The above parameters are given by GPS data collecting software and become from site calibration using pairs of know [local cs] and surveyed [wgs84] points.

Example input coordinates:

latitude = 54.3324735043758
longitude = 18.6188610214872

Expected result:

x = 30395.719
y = 28273.712

Whether is this possible only with those parameters provided or I need something else?
Can it be done using Proj4 or pyproj?
If yes, what would be apropriate command?
If not, can you hint me how to do it by hand, like links to formulas etc, please?

Thanks in advance.

[edit] Adding some background to my question.

In the case if my question seems absurd...
As I mentioned above the parameters become from site calibration or GPS alignmnet using pairs of know and surveyed points. I was told that after alignment is done, what I have to do to obtain x,y coordinates for any next surveyed wgs84 point, is to apply stereographic projection of WGS84 ellipsoid on a plane where:

B0, L0 are geodetic coordinates of the tangent point
X0, Y0 are grid coordinates of the tangent point
scale [k] is the scale factor

and with the least square method parameters of the transition to the local cs are determined.

Sorry if anything here is unclear but I am not quite familiar with subject.
Any help would be appreciated.

[edit]Adding additional sample points.

latitude longitude X Y
54.438657278 18.4602869845 41989.871 17718.469
54.4192483588 18.4806069996 39857.12 19081.45
54.4124384305 18.4716474826 39087.23 18515.701
54.3891511496 18.4799499633 36506.74 19108.41
54.3410205023 18.5443421778 31239.26 23406.7
54.3498711013 18.647572105 32374.58 30096.13
54.3515138914 18.6760425965 32600.628 31942.594
54.3800402814 18.6055277359 35669.507 27287.217
54.3971889473 18.6178833183 37596.017 28046.484
54.4041298674 18.6607429298 38432.47 30811.43
54.4070050071 18.6375394617 38717.65 29297.69
54.3609834076 18.474278312 33364.349 18804.941

[edit]Final word. Lests start from the beginning.
It turned out that my question was not quite right, in its original form.
So I decided to post full context and all steps I make in order to achieve the ultimate goal which is: getting local grid coordinates.

1.Site calibration points pairs.

latitude longitude x y
54.4353877827 18.451412164 41614.22 17150.17
54.3625019774 18.4469405743 33496.96 17024.81
54.3097484306 18.5670103377 27791.39 24957.15
54.3084348677 18.6299997771 27737.21 29059.27
54.402430173 18.6683595667 38254.91 31310.35
54.4230519006 18.6017101505 40450.41 26931.58
54.438657278 18.4602869845 41989.871 17718.469
54.4192483588 18.4806069996 39857.12 19081.45
54.4124384305 18.4716474826 39087.23 18515.701
54.3891511496 18.4799499633 36506.74 19108.41
54.3410205023 18.5443421778 31239.26 23406.7
54.3498711013 18.647572105 32374.58 30096.13
54.3515138914 18.6760425965 32600.628 31942.594
54.3800402814 18.6055277359 35669.507 27287.217
54.3971889473 18.6178833183 37596.017 28046.484
54.4041298674 18.6607429298 38432.47 30811.43
54.4070050071 18.6375394617 38717.65 29297.69
54.351702123 18.6831419293 32632.47 32403.62
54.4427340176 18.4391134326 42415.52 16335.76
54.4427368643 18.4592502382 42442.44 17641.9
54.3609834076 18.474278312 33364.349 18804.941
54.3714455799 18.7279691735 34899.27 35264.2
54.2911268152 18.6361668642 25820.32 29504.89

2.Convert above surveyed latitude and longitude to stereographic plane coordinates, taking coordinates of first calibration point as central latitude and longitude.
I don't know how important the scale factor [k_0] is so decided to force 1.0.
I also experimented with UTM instead of stereographic plane but result was not satisfactory with residuals = ~0.1m because we are talking here about obtaining survey grade coordinates not gis accuracy.

import pyproj

#input coordinates    
p1 = pyproj.Proj(proj="latlong",datum="WGS84", ellps="WGS84")

#output result
p2 = pyproj.Proj(proj="stere", lat_0=54.4353877827032, lon_0=18.4514121640352, k_0=1.0, x_0=0, y_0=0, ellps="WGS84")

stereographic_coords= pyproj.transform(p2, p1, latitude, longitude)
print stereographic_coords[1], stereographic_coords[0] #keep N(=x), E(=y) order

The result:

N E
3.50449394074e-07 2.28385549832e-06
-8113.07087572 -290.673432675
-13978.9604246 7524.00783763
-14116.6372391 11624.2036175
-3646.95639951 14088.88636
-1362.75920089 9755.69391324
363.973973779 575.837791136
-1796.13420238 1895.18404979
-2554.36871405 1313.79387579
-5146.3420395 1853.8865397
-10500.2307161 6044.01338437
-9501.32640299 12755.170692
-9312.95200167 14605.8597464
-6149.937535 10013.9390815
-4239.26761866 10812.2644791
-3459.25317782 13593.690851
-3143.43494621 12086.0443555
-9290.50857805 15067.4007146
817.802325033 -797.917280294
818.077929071 508.518662863
-8281.8734697 1486.45713731
-7082.36334017 17973.5167497
-16042.1577585 12030.6541087

3.Similarity transformation.
We can perform now similarity transformation to obtain local coordinates for any given point. One thing to remember is we firstly need to convert given latitude, longitude to stereographic coordinates.

transformation.py

import math
import numpy as np


class GeometricTransform(object):
"""Perform geometric transformations on a set of coordinates.

"""
    def __call__(self, coords):
       raise NotImplementedError()

    def inverse(self, coords):
        raise NotImplementedError()

    def __add__(self, other):
        raise NotImplementedError()


class ProjectiveTransform(GeometricTransform):
    _coeffs = range(8)

    def __init__(self, matrix=None):
        if matrix is None:
        # default to an identity transform
            matrix = np.eye(3)
        if matrix.shape != (3, 3):
            raise ValueError("invalid shape of transformation matrix")
        self._matrix = matrix

    @property
    def _inv_matrix(self):
        return np.linalg.inv(self._matrix)

    def _apply_mat(self, coords, matrix):
        coords = np.array(coords, copy=False, ndmin=2)

        x, y = np.transpose(coords)
        src = np.vstack((x, y, np.ones_like(x)))
        dst = np.dot(src.transpose(), matrix.transpose())

    # rescale to homogeneous coordinates
        dst[:, 0] /= dst[:, 2]
        dst[:, 1] /= dst[:, 2]

        return dst[:, :2]

    def __call__(self, coords):
        return self._apply_mat(coords, self._matrix)

    def inverse(self, coords):
        return self._apply_mat(coords, self._inv_matrix)

    def estimate(self, src, dst):
        xs = src[:, 0]
        ys = src[:, 1]
        xd = dst[:, 0]
        yd = dst[:, 1]
        rows = src.shape[0]

    # params: a0, a1, a2, b0, b1, b2, c0, c1
        A = np.zeros((rows * 2, 9))
        A[:rows, 0] = xs
        A[:rows, 1] = ys
        A[:rows, 2] = 1
        A[:rows, 6] = - xd * xs
        A[:rows, 7] = - xd * ys
        A[rows:, 3] = xs
        A[rows:, 4] = ys
        A[rows:, 5] = 1
        A[rows:, 6] = - yd * xs
        A[rows:, 7] = - yd * ys
        A[:rows, 8] = xd
        A[rows:, 8] = yd

    # Select relevant columns, depending on params
        A = A[:, self._coeffs + [8]]

        _, _, V = np.linalg.svd(A)

        H = np.zeros((3, 3))
    # solution is right singular vector that corresponds to smallest
    # singular value
        H.flat[self._coeffs + [8]] = - V[-1, :-1] / V[-1, -1]
        H[2, 2] = 1

        self._matrix = H

    def __add__(self, other):
    """Combine this transformation with another.

    """
        if isinstance(other, ProjectiveTransform):
        # combination of the same types result in a transformation of this
        # type again, otherwise use general projective transformation
            if type(self) == type(other):
                tform = self.__class__
            else:
                tform = ProjectiveTransform
            return tform(other._matrix.dot(self._matrix))
        else:
             raise TypeError("Cannot combine transformations of differing "
                        "types.")





class SimilarityTransform(ProjectiveTransform):
"""2D similarity transformation of the form::

        X = a0*x - b0*y + a1 =
          = m*x*cos(rotation) + m*y*sin(rotation) + a1

        Y = b0*x + a0*y + b1 =
          = m*x*sin(rotation) + m*y*cos(rotation) + b1

where ``m`` is a zoom factor and the homogeneous transformation matrix is::

    [[a0  b0  a1]
     [b0  a0  b1]
     [0   0    1]]

Parameters
----------
matrix : (3, 3) array, optional
    Homogeneous transformation matrix.
scale : float, optional
    Scale factor.
rotation : float, optional
    Rotation angle in counter-clockwise direction as radians.
translation : (tx, ty) as array, list or tuple, optional
    x, y translation parameters.

"""

    def __init__(self, matrix=None, scale=None, rotation=None,
             translation=None):
        params = any(param is not None
                 for param in (scale, rotation, translation))

        if params and matrix is not None:
            raise ValueError("You cannot specify the transformation matrix and "
                         "the implicit parameters at the same time.")
        elif matrix is not None:
            if matrix.shape != (3, 3):
                raise ValueError("Invalid shape of transformation matrix.")
            self._matrix = matrix
        elif params:
            if scale is None:
                scale = 1
            if rotation is None:
                rotation = 0
            if translation is None:
                translation = (0, 0)

            self._matrix = np.array([
                [math.cos(rotation), - math.sin(rotation), 0],
                [math.sin(rotation),   math.cos(rotation), 0],
                [                 0,                    0, 1]
            ])
            self._matrix *= scale
            self._matrix[0:2, 2] = translation
        else:
        # default to an identity transform
            self._matrix = np.eye(3)

    def estimate(self, src, dst):
        xs = src[:, 0]
        ys = src[:, 1]
        xd = dst[:, 0]
        yd = dst[:, 1]
        rows = src.shape[0]

    # params: a0, a1, b0, b1
        A = np.zeros((rows * 2, 5))
        A[:rows, 0] = xs
        A[:rows, 2] = - ys
        A[:rows, 1] = 1
        A[rows:, 2] = xs
        A[rows:, 0] = ys
        A[rows:, 3] = 1
        A[:rows, 4] = xd
        A[rows:, 4] = yd

        _, _, V = np.linalg.svd(A)

    # solution is right singular vector that corresponds to smallest
    # singular value
        a0, a1, b0, b1 = - V[-1, :-1] / V[-1, -1]

        self._matrix = np.array([[a0, -b0, a1],
                             [b0,  a0, b1],
                             [ 0,   0,  1]])

    @property
    def scale(self):
        if math.cos(self.rotation) == 0:
        # sin(self.rotation) == 1
            scale = self._matrix[0, 1]
        else:
             scale = self._matrix[0, 0] / math.cos(self.rotation)
        return scale

    @property
    def rotation(self):

        return math.atan2(self._matrix[1, 0], self._matrix[1, 1])

    @property
    def translation(self):
        return self._matrix[0:2, 2]




TRANSFORMS = {
    'similarity': SimilarityTransform,
}
HOMOGRAPHY_TRANSFORMS = (
    SimilarityTransform,
    ProjectiveTransform
) 


def estimate_transform(ttype, src, dst, **kwargs):
    ttype = ttype.lower()
    if ttype not in TRANSFORMS:
        raise ValueError('the transformation type \'%s\' is not'
                     'implemented' % ttype)

    tform = TRANSFORMS[ttype]()
    tform.estimate(src, dst, **kwargs)

    return tform


def matrix_transform(coords, matrix):
    return ProjectiveTransform(matrix)(coords)

And here's my calculation:

import numpy as np
import transformations as tf #importing above script

##############################
####CALIBRATION POINTS########
##############################

#primary system coordinates obtained instep 2 from surveyed latitude,longitude.

ps1=( 3.50449394074e-07 , 2.28385549832e-06 )
ps2=( -8113.07087572 , -290.673432675 )
ps3=( -13978.9604246 , 7524.00783763 )
ps4=( -14116.6372391 , 11624.2036175 )
ps5=( -3646.95639951 , 14088.88636 )
ps6=( -1362.75920089 , 9755.69391324 )
ps7=( 363.973973779 , 575.837791136 )
ps8=( -1796.13420238 , 1895.18404979 )
ps9=( -2554.36871405 , 1313.79387579 )
ps10=( -5146.3420395 , 1853.8865397 )
ps11=( -10500.2307161 , 6044.01338437 )
ps12=( -9501.32640299 , 12755.170692 )
ps13=( -9312.95200167 , 14605.8597464 )
ps14=( -6149.937535 , 10013.9390815 )
ps15=( -4239.26761866 , 10812.2644791 )
ps16=( -3459.25317782 , 13593.690851 )
ps17=( -3143.43494621 , 12086.0443555 )
ps18=( -9290.50857805 , 15067.4007146 )
ps19=( 817.802325033 , -797.917280294 )
ps20=( 818.077929071 , 508.518662863 )
ps21=( -8281.8734697 , 1486.45713731 )
ps22=( -7082.36334017 , 17973.5167497 )
ps23=( -16042.1577585 , 12030.6541087 )

#secondary system coordinates/known local coordinates
ss1 = ( 41614.22 , 17150.17 )
ss2 = ( 33496.96 , 17024.81 )
ss3 = ( 27791.39 , 24957.15 )
ss4 = ( 27737.21 , 29059.27 )
ss5 = ( 38254.91 , 31310.35 )
ss6 = ( 40450.41 , 26931.58 )
ss7 = ( 41989.871 , 17718.469 )
ss8 = ( 39857.12 , 19081.45 )
ss9 = ( 39087.23 , 18515.701 )
ss10 = ( 36506.74 , 19108.41 )
ss11 = ( 31239.26 , 23406.7 )
ss12 = ( 32374.58 , 30096.13 )
ss13 = ( 32600.628 , 31942.594 )
ss14 = ( 35669.507 , 27287.217 )
ss15 = ( 37596.017 , 28046.484 )
ss16 = ( 38432.47 , 30811.43 )
ss17 = ( 38717.65 , 29297.69 )
ss18 = ( 32632.47 , 32403.62 )
ss19 = ( 42415.52 , 16335.76 )
ss20 = ( 42442.44 , 17641.9 )
ss21 = ( 33364.349 , 18804.941 )
ss22 = ( 34899.27 , 35264.2 )
ss23 = ( 25820.32 , 29504.89 )

############################
############################
############################

from_pt = [ps1, ps2, ps3, ps4, ps5, ps6, ps7, ps8, ps9, ps10, ps11, ps12, ps13, ps14,  ps15, ps16, ps17, ps18, ps19, ps20, ps21, ps22, ps23]
to_pt =   [ss1, ss2, ss3, ss4, ss5, ss6, ss7, ss8, ss9, ss10, ss11, ss12, ss13, ss14, ss15, ss16, ss17, ss18, ss19, ss20, ss21, ss22, ss23]

fr = np.array(from_pt)
to = np.array(to_pt)
tform = tf.estimate_transform('similarity', fr, to)

#list of coordinates to transform, first converted to sterographic
l = [
[10892.8449025, -11442.6688292],
[10892.2806078, -11442.7145534],
[10880.2426997, -11440.0353547],
[10869.0233644, -11437.2787524],
[10868.4894383, -11437.3154292],
[10879.4637253, -11439.8965462],
[10879.194897, -11440.0837157],
[11661.9477535, -12607.9800965],
[11620.0600416, -12597.7915158],
[11600.121001, -13934.58708]]

for line in l:
    x = line[1]
    y = line[0]
    res= tform(np.array([x,y]))
    print res

Final result:

my calculation result                   known control points coordinates
[[ 30395.74138597  28273.70706175]]     # 30395.716 28273.716
[[ 30395.68418212  28273.14381546]]     # 30395.659 28273.153
[[ 30398.1177287   28261.05386387]]     # 30398.092 28261.063
[[ 30400.64532909  28249.78073875]]     # 30400.620 28249.790
[[ 30400.59778907  28249.24767053]]     # 30400.572 28249.257
[[ 30398.24064822  28260.27222544]]     # 30398.215 28260.281
[[ 30398.04804427  28260.0072639 ]]     # 30398.022 28260.016
[[ 29246.33177298  29066.37582526]]     # 29246.308 29066.384
[[ 29255.66538993  29024.28939019]]     # 29255.641 29024.297
[[ 27918.74210572  29031.57195305]]     # 27918.718 29031.581

I think, residuals are acceptable.

Answer 1

For such problems, I try to visualize the points in QGIS to see where they are placed.

From your parameters, I created a custom CRS with the definition:

+proj=sterea +lat_0=54.4353877827032 +lon_0=18.4514121640352 +k=0.999790760649094 +x_0=41614.2107651061 +y_0=17150.1692507701 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

and created points in WGS84 CRS and the Stereo projection. Unfortunately, they do not align. Transforming the WGS84 coordinates to your Stereo projection I get:

    X: 52504.84
    Y: 5709.89

Visualizing your dataset, I get the following picture:

So what went wrong? You simply swapped X and Y. "Normally", we think of Easting as X and Northing as Y. Your parameters are right, but you have to swap X0 and Y0:

+proj=sterea +lat_0=54.4353877827032 +lon_0=18.4514121640352 +k=0.999790760649094 +x_0=17150.1692507701 +y_0=41614.2107651061 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

I added the points once in WGS84 and once in the Stereo projection. Looking closer at the picture, you see an offset between some points of around 100 metres. This is the point where usually the rotations from the +towgs84 parameters come into effect. Or the central meridian was not the right one.

So now we have a proj definition which you can use with GDAL cs2cs to convert the point coordinates using a batch with the command line:

cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=sterea +lat_0=54.4353877827032 +lon_0=18.4514121640352 +k=0.999790760649094 +x_0=17150.1692507701 +y_0=41614.2107651061 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs in.txt >>out.txt

for as many points as you want.

Doing some linear regression on my own, I found this projection string best fitting for your sample data:

+proj=tmerc +lat_0=54 +lon_0=17.02 +k=1 +x_0=-75731 +y_0=-7791 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

The cs2cs command line now is:

cs2cs +init=epsg:4326 +to +proj=tmerc +lat_0=54 +lon_0=17.02 +k=1 +x_0=-75731 +y_0=-7791 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs in.txt >out.txt

Don't bother that it is not a stereographic projection anymore: For the size of your area, the difference between projection methods is much smaller than finding the right central meridian.

For the linear regression, I used a simple LibreOffice Calc sheet:

Imported your sample data, and made an input file for cs2cs from the first two columns with easting first.

As first try, I chose lat_0 and lon_0 rounded to full degrees of the sample data setting x0 and y0 to zero, ran the cs2cs command and filled the output in column E to G. The next columns compute the difference between target and cs2cs values and the absolute offset. Below, I computed mean and standard deviation for deltaX deltaY and radius.

After that, I played a bit with the cs2cs parameters, and found that lat_0, x0 and y0 have no effect on the std dev value. So I changed lon_0 until the std dev was less than 1 meter. I copied the std dev row for each run as protocol.

The mean values with inverted sign then give the values for x0 and y0 for the final step where delta X and deltaY is less than 2 metres for all points.

You can do the same with different (but similar) projections, values for lat_0, k, ellipsoid or +towgs84. The sample points should cover your area of interest, as extrapolation will give worse results.

---

Final Edit

Taking the tangent point from polish UKLAD 65 Zone 3 as mentioned in the comments, I got best results with the following command line:

cs2cs +init=epsg:4326 +to +proj=sterea +lat_0=53.58333333333334 +lon_0=17.00833333333333 +k=0.9998 +x_0=-76625.57 +y_0=-54048.732 +ellps=krass +towgs84=33.4,-146.6,176,0,0,0,0 +units=m +no_defs in.txt >out.txt

Offset is less than 1 meter for all sample points, and resulting standard dev 0.32m. You could even try to use all 7 Helmert parameters to get it better.

Answer 2

The "localization" performed by GNSS controllers use a few different methods to define a least-squares fitting coordinate system.

In nearly all cases, a temporary Transverse Mercator projection on a WGS84 ellipsoid is used. The central meridian of the TM proj is the center-of-gravity of the input control points used in the adjustment. Then a rotation, translation and scaling are performed about that point. These are just the horizontal values, as vertical values are handled seperately.

So you'll end up with 6 parameters:

1. Origin Lat
2. Origin Lon
3. Rotation about Origin
4. Translation (x) from Origin
5. Translation (y) from Origin
6. Scale

I'm guessing one of your mystery values (q) is the rotation about the origin...

Hope this helps.