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.
Best Answer
(My original answer forgot about your position and radius being in different units!)
Assuming accuracy is unimportant, then if your centre is (X, Y) and your radius is R, the bounding box corners are simply (X-R, Y-R) and (X+R, Y+R).Considering the different units:
Assuming accuracy is unimportant, then if
(X, Y)
, whereY
is latitude,X
is longitude, in degreesr
, andR
, in the same units asr
Then the bounding box corners are simply
(X - dX, Y - dY)
and(X + dX, Y + dY)
wheredY = 360 * r / R
, the "search radius", as difference in latitude,dX = dY * cos (rad (Y))
, the "search radius", as difference in longitude.