[Proj] Add simple transformations to proj

[Jose Gonçalves]

Hello

I had a similar problem with some old scanned maps from Portugal that were
in a local topographic coordinate system. I also wanted to have a PROJ.4
string for that system in order to load the maps in a GIS software together
with other data sets. I did it simulating a LCC projection (Lambert
Conformal Conic).

First there is an assumption we can do about the affine transformation. It
is likely that it is nearly a conformal transformation (a=e and d=-b) and
with a scale factor very close to 1 (sqrt(a*a+b*b)=1). In this case the
essential thing to deal with is the rotation of the axes.

In the LCC projection you can choose one standard parallel in your region
so that the linear deformation is always very close to 1. Then you can
choose a central meridian far enough, so that the cartographic north (the Y
axis of the projection) becomes parallel to the Y axis of your local
coordinate system. You can do that by trial and error or by analysing the
meridian convergence of the LCC projection. Then you can add a false
easting/northing so that the projection and your local system coincide.

Suppose for example that point (4.8E, 52.2N) is the origin of your local
coordinate system and point (4.9E, 52.4N) is on the Y axis.
Consider 52.3N as the standard parallel. The following PROJ string converts
these two points in a way that they have the same easting:

   +ellps=bessel +proj=lcc +lon_0=-16.68875 +lat_0=52.3 +lat_1=52.3
+lat_2=52.3 +no_defs

The projected coordinates of the first point are (1447575.65 m, 205246.75
m). If you add symmetric values as false easting/northing

   +x_0=-1447575.65 +y_0=-205246.75

the PROJ string does what we want. You can also play with k_0 if some scale
adjustment is still needed.

This simulation of a projection is obviously only valid for rather small
areas. It worked well for me, in Quantum GIS, with the "on the fly" option.


Regards

Jose Gonçalves



2013/1/2 Jan Hartmann <j.l.h.hartmann at uva.nl>

>  Hi Frank and all, happy new year!
>
> Is the following possible with Proj? I am working with historical maps
> that are rotated with respect to a well known epsg-projection, see for
> example
>
> http://mapserver.sara.nl/a1000/bladgrenzen.html
>
> The blue squares are the sheets of the map 1:1000 of Amsterdam after 1940,
> in epsg:28992, the red squares are the same sheets before 1940 in a local
> coordinate system. I know the conversion coefficients for the two
> coordinate systems: a simple six-value transformation of the form:
>
> xnew = a*xold + b*yold + c
> ynew = d*xold + e*yold + f
>
> I would like to access the second map via a proj-string. For that, I would
> have to expand proj with a procedure that takes the six coefficients (a-f)
> and generates the transformed coordinates after the regular transformation.
>
> Is this possible? Can anyone give me some pointers where to look in the
> proj-code? Conceptually, it looks simple to me, but I don't know where to
> start.
>
> Jan
>
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20130103/5b89d7f6/attachment.html>