How can I obtain a natural projection of a raster map from southern Poland that displays equal north-south and east-west distances?

Question

I am attempting to make a digitised version of an 1848 cadastral map of Siemiechów, Poland (Lat 49.853, Lon 20.902). I have 44 map images that I aligned and assembled into a large mosaic. I used the QGIS raster georeferencer and a selection of latitude/longitude data for road intersections to prepare a georeferenced TIF file. I used various CRSs (see below). I am using QGIS 3.22.14 desktop.

No matter which CRS I choose, the georeferenced on-screen image is squashed, that is, the north-south dimensions/distances appear much smaller than east-west. This is a problem with how the raster is projected. I have tried:

EPSG 2178 – ETRF2000-PL / CS2000/21

EPSG 3035 – ETRS89-Extended /LAEA Europe

EPSG 3857 – WGS84 Pseudo Mercator

EPSG 4326 – WGS84

ESRI 102176 – ETRS_1989_UWPP_2000_PAS_7

I checked my work by snipping map images from Google and OpenStreetMap and georeferencing them. The same substantially squashed georeferenced raster projections were obtained from maps from both sources.

I understand that every projection distorts lengths. But the distortions I'm seeing are too extreme. Is there a projected CRS that would (roughly) replicate the X-Y dimensions of the original map? My mosaic raster is much larger than 1 Gb, so I show the same issue with an OpenStreetMap raster from the same region as a proxy for my problem (below).

Do I need to use a custom coordinate reference system to obtain a naturally projected georeferenced map image/raster that displays (nearly) equal north-south distances and east-west distances? Or is there an existing CRS that will do this? Is there an existing solution in QGIS that can achieve this?

Here is an example showing an original raster from OpenStreetMap and the georeferenced TIF.

I used the following reference points:

Point Lon Lat (decimal degrees)

1 20.875906 49.856044

2 20.895682 49.843258

3 20.936216 49.844372

4 20.892439 49.87169

And here is the georeferenced raster using EPSG:4326 - WGS84:

Answer 1

The answer is to use ETRS89-extended / LAEA Europe -OR- ETRF2000-PL / CS2000/21 to georeference and display the raster mosaic. I came to this result by extracting the relevant map image from OpenStreetMap and using the QuickOSM Plugin to display the boundaries in QGIS. I then used the QGIS View > Decorations > Grid command to show that vertical and horizontal distances were equal. From there I viewed the vertices of the town boundaries. Not only did this give me confidence that I was displaying a natural map, it also gave me lots of points to georeference my raster.

Earlier I came up with a simple work-around, admittedly not advised. I georeferenced a screen-shot of an OpenStreetMap raster of the same area using EPSG:3857. The reference points I used were near the North-South and East-West borders of the area of interest. When the dimensions of the georeferenced TIF raster were compared to those of the original raster, I found that the latitude dimension was reduced by a factor of 0.646696035. I used this to rescale the reference latitudes using the formula Lmin +(L-Lmin)/SF, where Lmin=minimum latitude, L=latitude value to be scaled, SF=scaling factor. When I georeferenced the raster again using scaled latitudes, the latitude dimension was reduced by less than 0.2%, which is undetectable by eye. In other words, the projection of the raster was restored to virtually the original perspective that was displayed by OpenStreetMap. I could use this scaling on the raster I want to convert to a vector map. It turns out that the scaling factor is Cos(Latitude).