guessproj
is a Python script that calculates unknown parameters
of cartographic projection or coordinate system from coordinates
of identical points in some known coordinate system and in unknown one.
You should know projection type, though.
The script can also determine parameters of transformation between two datums.
The script uses pyproj
and scipy
internally.
The method of least squares is used, so the more points you have,
the better accuracy will be achieved.
Python 2.6+ and 3.3+ are supported. The script is written in pure Python itself but depends on some packages that are not.
The best way to install guessproj
is using
pip:
pip install guessproj
Be aware that guessproj
has some binary dependencies that you need
to install before trying to install guessproj
.
These are GDAL and PROJ.4.
If you are a GIS specialist you probably already have these libraries.
You also need Python bindings for them (GDAL bindings and pyproj)
as well as NumPy and SciPy packages.
The pip
tool will try to install these packages automatically
but on most systems you'll need to install them in a platform-specific way.
Note that GDAL is an optional dependency but may become required in future
versions of guessproj
.
Instead of using pip
, you can download the source archive,
unpack it and run
python setup.py install
in the unpacked directory. You need setuptools to do this.
Also, you can use guessproj
without installation. Just download
the file guessproj.py and run it like any Python script:
python guessproj.py --help
You still need dependencies to be installed, of course.
The input data for the script is a text file that contains lines of space separated values, each line representing a point.
Input data format for 2D points is as follows:
x1 y1 x2 y2 point name
Input data format for 3D points is as follows:
x1 y1 z1 x2 y2 z2 point name
Point name may contain spaces but should not start with a number.
If the line starts with #
, it will be ignored.
The coordinates in known coordinate system go first, followed by coordinates in unknown one. 3D points are mostly useful if there exists a datum transformation between coordinate systems.
Obviously, the number of point coordinate values must be not less than the number of unknown parameters to determine.
Both decimal format and degrees, minutes and seconds format
(XXXdXX'XX.XXX"
) are supported for coordinate values.
You can use comma as well as period to separate fractional part.
Command line syntax used to run the script is similar to that of cs2cs
utility which comes with PROJ.4
library. You should specify parameters
of known coordinate system in PROJ.4
format (referred as projstring here),
known and unknown parameters of the unknown system, and a path to input file.
Exact parameters of various coordinate systems you can find
in PROJ.4 documentation.
The unknown parameters are specified among all others in the projstring,
the only difference is using ~
symbol instead of =
. The numeric value
that follows the ~
symbol is an initial approximation of the parameter value.
For +towgs84
parameter, you can specify ~
before any (or all) of comma-separated
values. The combination =~
is the same as ~
.
Example:
guessproj +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to \ +proj=tmerc +ellps=krass +lat_0=0 +lon_0~44 +x_0=300000 +y_0~-4.7e6 \ +towgs84=23.57,-140.95,-79.8,0,0.35,0.79,-0.22 points.txt
All that goes before +to
argument is a projstring for the known system
(if omitted, WGS84 longitude/latitude is used by default). All the rest
parameters starting with +
are the projstring for the unknown system,
where initial approximations of the unknown parameters are marked with ~
.
In this example, parameters +lon_0
and +y_0
are unknown. The last argument
is a name of input text file containing point coordinates.
The script can evaluate numeric parameters only, so you should specify
at least +proj
and +ellps
. It's worth mentioning, also, that some
parameters have the same (or nearly the same) effect, so it's a bad idea,
for example, to specify both +lat_0
and +y_0
as unknown
for Transverse Mercator projection. The system of equations will be
ill-determined in that case.
Any additional program options (which are not part of projstring syntax)
start with -
or --
.
Option -h
or --help
prints a short command line reference and exits
the program.
Option --encoding=ENCODING_NAME
specifies the encoding of input file
(--encoding=utf-8
by default).
Option --proj
or --proj4
forces output of resulting projstring only,
suppressing table of residual errors.
Option --wkt
forces output of projection parameters in OGC WKT format
(GDAL Python bindings required).
Option --esri
forces output of projection parameters in Esri WKT format
(experimental, GDAL Python bindings required).
Option --pretty
forces pretty WKT formatting when used with --wkt
or --esri
.
The default output of the program is a projstring in which approximated values of parameters are replaced with the exact values found by the script, and a list of residual errors for each point. Other forms of output can be specified using program options.
If GDAL bindings are installed, the projstring will be formatted so as to be represented in a normalized form.
To run unit tests with Python 2.7 or 3.3+, execute in source directory:
python -m unittest discover test
In Python 2.6, you should install unittest2 package and use:
PYTHONPATH=. unit2 discover test
You can also run scripts from test/
directory directly.