Chris Crook | 5 Jan 09:01 2015

Time dependent coordinate transformations

I'm wanting open a discussion around how to handle time dependent coordinate transformations.  These are
becoming far more relevant to GIS data sets and positioning than they have historically been.

This is a growing issue for which I'd like to find the right forums and process to ultimately have have widely
supported in coordinate transformation implementations, so please suggest or forward other
appropriate forums.

I'm posting initially to the PROJ list as PROJ underlies much of the coordinate system transformations in
the open source GIS stack, and does include datum transformations as well as projections.  However this is
clearly a wider issue .. as I suggest below I think this affects the definition of spatial features
generally - so I'd welcome suggestions for more appropriate or additional forums.

The background

Coordinate systems are based on geodetic datums.  These are generally either international datums, such
as WGS84 or ITRF realisation, or national datums defined by a country's geodetic authority.

Within international datums apparently fixed features (such as buildings) are moving due to plate
tectonics at rates of typically 5 cm/year.  This means that the WGS84 coordinate of a "fixed" point in 2000
would now miss that point by 75cm.

National datums are generally referenced to fixed marks within the country and provide a stable
coordinate system, such that the coordinates of fixed marks are constant.

Obviously the conversion between international datums and national datums cannot be done without taking
account of when the coordinates are applicable.

Typically GIS systems are in terms of national coordinate systems.  Generally they contain spatial
definitions in terms of the national datum, with no timestamp, and assume that this provides a consistent
(Continue reading)

Charles Karney | 20 Dec 15:33 2014

Rhumb lines and Mercator on a triaxial ellipsoid

It is well known that a rhumb line arrives at a pole in a finite
distance after encircling the pole infinitely many times.  Craig Rollins
recently asked me what heading a rhumb line has *after* passing through
the pole.

One way of answering this is to consider a rhumb line on a triaxial
ellipsoid and to take the limit as the two large axes approach one
another.  The (somewhat surprising) result is that the heading of the
rhumb line is reversed.  E.g., if the initial heading is NE, the heading
after passing through the pole is SW.

This got me to thinking about the Mercator projection on a triaxial
ellipsoid.  This was given by Jacobi in 1843, see section 28 of

The integrals that Jacobi gives can be written in terms of elliptic
integrals, see

Finally, Jacobi has an interesting take on Gauss' work on conformal
projections (excerpted from Balagangadharan's translation):

   "Among the different ways of representing a curved surface on a plane,
   as is necessary for a map, one prefers, above all, the method of
   projection in which infinitely small elements remain similar.  In the
   preceding century Lambert had been concerned with various aspects of
   this projection, of which one can learn in detail from his
(Continue reading)

Frank Willett | 18 Dec 23:11 2014

Release Estimate?

Anyone have a guess as to when the new release will come? I'm waiting to
update projects which use proj.

Frank Willett
Proj mailing list
Proj <at>

Margherita Di Leo | 17 Dec 14:34 2014

convert latlon to UTM without knowing the zone beforehand


I need to convert any pairs of geographic coordinates (latlon) into UTM. This needs to run automatically so that I need an automatic tool to tell for any pairs which UTM zone it belongs to before converting. Of course I can encode in my script all possible cases, but before reinventing the wheel i wanted to ask if there is such tool already within proj.4 library.

Thank you in advance

Best regards,

Dr. Margherita DI LEO    
Scientific / technical project officer

European Commission - DG JRC 
Institute for Environment and Sustainability (IES)
Via Fermi, 2749
I-21027 Ispra (VA) - Italy - TP 261
Tel. +39 0332 78 3600   

Disclaimer: The views expressed are purely those of the writer and may not in any circumstance be regarded as stating an official position of the European Commission.
Proj mailing list
Proj <at>
Charles Karney | 16 Dec 14:14 2014

The area of a rhumb line polygon

The calculation of the area of a geodesic polygon has been included in
GeographicLib since v1.4 (2010-09).

In the latest version v1.39, I've also included the equivalent
calculation when the edges of the polygon are rhumb lines.  The
formulation is given in

and the online planimeter

implements this area calculation.  I've also worked out the area when
the edges are great ellipses.  For the formulation, see

An implementation is included in the Matlab package

Comments are welcome.

Proj mailing list
Proj <at>

SriRam Prasad Bhasker | 14 Nov 07:24 2014

Build GDAL with NADSHIFT's


I have a below .prj file


GEOGCS["GCS_NAD_1983_2011",DATUM["D_NAD_1983_2011",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433],METADATA["USA - CONUS and Alaska; PRVI",167.65,14.93,-63.89,74.71,0.0,0.0174532925199433,0.0,1511]]


Can someone guide me how to project a Shape file with the below prj File.

Testepsg is giving a failure for the above prj file.


Let me know if anyone has used this earlier… To use the Conus and Alaska Parameters

"Disclaimer: This message and any attachments contain information that may be RMSI Private Limited. confidential and/or privileged. If you are not the intended recipient (or authorized to receive for the intended recipient), and have received this message in error, any use, disclosure or distribution is strictly prohibited. If you have received this message in error, please notify the sender immediately by replying to the e-mail and permanently deleting the message from your computer and/or storage system"
Proj mailing list
Proj <at>
Jeff McKenna | 13 Nov 13:14 2014

Re: PROJ 4.9.0RC1

I may have missed the public message, but in case others did as well:
there is an RC2 also available for testing:


Jeff McKenna
MapServer Consulting and Training Services

On 2014-09-13 8:33 PM, Frank Warmerdam wrote:
> Folks,
> I have discovered that PROJ.4 was left with a hanging beta
> again some time ago without my following up on the release.
> Grr, I'm lame!  Anyways, Howard Butler and I have made a
> pass through open bugs and fixed some.  Howard has also
> introduced CMake support for PROJ.4 (though the automake
> way is also still supported).
> Source is at:
> I would appreciate feedback before finalizing the release in the coming
> week.
> Best regards,
> Frank

Proj mailing list
Proj <at>

Shimi Chen | 5 Nov 10:47 2014

cs2cs produces partial output when called from matlab


I realize that MATLAB is proprietary software and thus not very reproducible.
However I would appreciate if someone could think if he has an idea why this could happen.
I will explain the issue with a simple example.
I made a file 0.tmp with the contents:

35.049 31.728
Running this command:
cs2cs +init=epsg:4326 +to +init=epsg:2039 <0.tmp >1.tmp Produces this correct 1.tmp file:
204718.31 626160.79 -20.79 However, when calling the exact same command from the matlab shell:
system(['cs2cs +init=epsg:4326 +to +init=epsg:2039 <0.tmp >1.tmp']) The resulting file 1.tmp only has these contents:
-20.787 So, clearly, cs2cs is getting all of the input, since otherwise it could not get the -20.79 right. The question is why would it not also write 204718.31    626160.79  to the file?
What weird thing could MATLAB be doing with the input that would make cs2cs produce this partial output?

Note: do not try to reproduce in Octave. It works fine in Octave.

Proj mailing list
Proj <at>
Zhang Qun | 25 Oct 08:46 2014

the x, and y meaning in pj_transform

Dear all

I have got a small question that I can not find the confirmed answer on the web. I am now using pj_transform to transform coordinates from svy21 to latlon using the following code snippet:

pj_svy21 = pj_init_plus("+proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs")
pj_latlong = pj_init_plus("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

double svy_east = 12056.915;
double svy_north = 37621.665;
double lat_radian = svy_east;
double lon_radian = svy_north;
projResult = pj_transform(pj_svy21, pj_latlong, 1, 0, &lat_radian, &lon_radian, NULL );

But the result (lat_radian, lon_radian) I got is reversed, i.e. the result given by
pj_transform(pj_svy21, pj_latlong, 1, 0, &x, &y, NULL )

is x corresponds to longitude and y corresponds to latitude. It is not the other way around as I did.

Is this the way it worked or I just used this function wrongly? Thanks.

Best regards,
Zhang Qun
Proj mailing list
Proj <at>
Artur Bercik | 24 Oct 08:26 2014

(no subject)

Dear Proj list members:

Being a very new in this field, could someone answer me about the following problem:

I have a pair of Geodetic Latitude and Geodetic Longitude (50.89524, 117.69043).
How can I convert it into Latitude and Longitude with WGS84 Ellepsoid?

Thanks in the advance.

Proj mailing list
Proj <at>
Charles Clercq | 23 Oct 10:13 2014

Convertion Lat/Long WGS84 -> UTM-WGS84 problem?

Dear all,

I have written a simple basic code inspired from the doc to convert data in lat/long WGS84 in meters (UTM/WGS84)
My problem is that when I tested it with lat: 51.0242 / long: 3.96338 and obtained  x: 7627600.282316 / y: 742130.588781 which I compared with online converter which gave me : 5652958 / 567564

I would like to understand the error. I that because proj can't convert from WGS84 lat/long to UTM WGS84 ?

I tested with command line:
echo 51.0242 3.96338 | cs2cs.exe +proj=latlong +ellps=WGS84 +to +proj=utm +ellps=WGS84 +datum=WGS84 +units=m
and obtained the exact same result.

Any advice?



The code
#include <iostream>

#include <proj_api.h>

int _tmain(int argc, _TCHAR* argv[])
projPJ pj_merc, pj_latlong;
    double x, y;

    if (!(pj_merc = pj_init_plus("+proj=utm +ellps=WGS84 +datum=WGS84 +units=m")) )
    if (!(pj_latlong = pj_init_plus("+proj=latlong +ellps=WGS84")) )
std::cout << "Lat: ";
double latitude;
std::cin >> latitude;
std::cout << std::endl;
std::cout << "Long: ";
double longitude;
std::cin >> longitude;

    x = latitude*DEG_TO_RAD;
y = longitude*DEG_TO_RAD;

    int p = pj_transform(pj_latlong, pj_merc, 1, 1, &x, &y, NULL );

//std::cout << "The coordinate in meters are: " << x << "/" << y << std::endl;
printf("Lat/Long: %f/%f - in meters ->: %f/%f\n", latitude, longitude, x, y);

return 0;
Proj mailing list
Proj <at>