Michal Seidl | 19 Apr 10:54 2014

Custom nadgrid and .lla format again

I would like to confirm my suspecting that there is a little mistake in
README.GRD file which describes *.lla format. As I can see values in lla
files distributed with proj it seems to me it is something like this

<line_no>: <x_shift> <y_shift> <dx_shift> <dy_shift> *

So there is line number, x and y shift, and then there are only differences 

I just guess from values. Second and next pair values are much smaller than
first pair of values.
There can be also problem in units. Are difference in micro (milionth)
seconds of arc? diff= diff(deg)*3600E6?

Here is what is documented

Subsequent lines are:

<line_no>: <x_shift> <y_shift> *

The <line_no> is zero based, and will vary from 0 to <grid_size_y>-1.  The 
number of x/y shift pairs will match <grid_size_x>.  Grid lines can be
split over multiple physical text lines.  Use the colon to identify starts
of new grid lines.  The shift values are in millionths of a secondc of arc. 

For example, from MD.lla:

(Continue reading)

raul pra levis | 26 Mar 15:16 2014

Add custom projection

<!-- .hmmessage P { margin:0px; padding:0px } body.hmmessage { font-size: 12pt; font-family:Calibri } -->
I work for a company where we have our custom mercator projection.
I already developed code snippets to integrate such projection in proj4 4.8.0  (c++).
Is there the chance to collaborate with develop team to integrate it in proj4 main trunk?

Thanks and best regards,

Raul Pra Levis

Proj mailing list
Proj <at> lists.maptools.org
Phil Elson | 4 Mar 15:41 2014

Expanding a proj4 CRS definition

I'm curious to know whether it is possible for proj4 to return a complete representation of a projection definition. For instance, supposing I gave a proj4 string of:


It would be useful to get back a projection definition such as:

+proj=merc +lat_ts=0 ...

Similarly, though this is less important to me, is it possible to go from:



+proj=aea +lat_1=34 +lat_2=40.5 ...

Without just returning the string found in share/proj/epsg?

I guess what I'm really asking whether there is an internal representation of these values which can then be used to pump out a string suitable for passing through to proj4 later on?

Thanks in advance,


Proj mailing list
Proj <at> lists.maptools.org
Frederik Ramm | 1 Mar 23:05 2014

Equal Arc Projection, GeoTIFF, and GDAL/Proj4


   I've been pointed to this list after posting my question to the
German FOSSGIS-talk list (FOSSGIS is the German OSGeo chapter).

Among the things I do for a living is making raster images from
OpenStreetMap data, usually through the "Mapnik" renderer which has
Proj4 support built-in. Someone recently asked whether I could make them
an image in the "Equal Arc" projection in GeoTIFF format, and I naively
said "sure, give me the EPSG code and I'll work it out". I had never
heard the term "Equal Arc" and I somewhat assumed that the inquirer
probably meant one of the common "Equal Area" projections.

But I was wrong; what the client-in-spe is looking for is indeed the
"Equal Arc" projection that is defined here


and which, as far as I can see, has neither an EPSG code nor anything
remotely machine readable like a WKT or Proj4 projection description.
Even so, it seems to be widely used in aviation, and it seems to be an
integral part of so-called ADRG files. These ADRG files can even be
"experimentally" written by GDAL:


but of course these aren't GeoTIFFs then. I wonder:

1. Is it at all technically possible to have a GeoTIFF in Equal Arc
projection, given that I could not find any represenation of that
projection suitable for including in a GeoTIFF header?

2. If the answer to (1) is "yes", then can I create such a GeoTIFF with
GDAL, or possibly even by passing a clever Proj init string to the
Mapnik renderer?

3. If the answer to (2) is "no", then is this just an issue where
someone would have to spend some time adding support to Proj (with any
luck I could get the client-in-spe to fund this work), or are there
fundamental reasons barring us from doing that?

My internet research regarding "Equal Arc" was unsatisfying; the name
"FalconView" appeared a number of times in conjunction with this

Thank you for any insights you might have on this.



Frederik Ramm  ##  eMail frederik <at> remote.org  ##  N49°00'09" E008°23'33"
Proj mailing list
Proj <at> lists.maptools.org

AdamDynamic | 24 Feb 10:37 2014

Beginner Question: Latitude/Longitude to x, y conversion - understanding output


I've been scouring the internet for a couple of days now and getting
absolutely nowhere - I appreciate this is probably a simple question once
you know the answer but any help on this would be appreciated!

I have a series of latitude/longitude coordinates that i want to plot as xy
coordinates (i.e. a flat, rectangular map).

I'm using the PyProj Python library to perform the conversion, the problem
is I can't determine what the output of the function means or represents? 

e.g. Converting [51.5072,0.1275] (what Google tells me is the lat/long for
London in the UK) returns [5733755.276187301, 14193.246790158479] for a
Mercator projection - my question is '5733755.276187301 of what'? I
understand that the library converts to Proj.4 notation (?) but I can't seem
to find what the units are, or how I'd map the output onto a 900x1100

I'm using the following conversion formula which I found  here

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
+k=1.0 +units=m  +no_defs

Any assistance would be great as I haven't been able to get any information
from anywhere else. Please let me know if more information would be useful.



View this message in context: http://osgeo-org.1560.x6.nabble.com/Beginner-Question-Latitude-Longitude-to-x-y-conversion-understanding-output-tp5105455.html
Sent from the PROJ.4 mailing list archive at Nabble.com.
Proj mailing list
Proj <at> lists.maptools.org

support.mn | 18 Feb 22:49 2014

More Proj.4 features (was Re: Aunits)


basically I have nothing against "aunits" as long as the default
behavior remains the same .. except that it is a very small
change and we would like to have more improvements
at the same time to the Proj.4 definition language.

One of the most important is the formal syntax checking and
check against most obvious user input errors (typos etc) ..
for example normal lat / lon limits on Earth etc. If errors were
detected the separate routine would then spit out a warning
or error message. This would then warn the user in fore
hand that there is an obvious error in the input. Let say
for example having an O (letter O) instead of a 0 (number 0)
in a number.. Now the sad situation is that Proj.4 says
nothing but gives crazy results and the user is not aware
of what is the problem and the crazy definition might
end up in some library and used with strange results.

The checking could be a separated subprogram that just
does the checking and nothing more. The main program
could just omit that with such data that is already checked
once to keep it simple and fast. So the check routine would
be a separated part and would not affect any of the current
library operation and programs that use it.

It would be most convenient to have that check inside
the library since it also defines the interface language
and format. Now everybody have to write it separated
if they want to have that.. situation similar to "aunits".

The syntax checker and analyzer could be also a
some kind of a preprocessor that could also handle
the "aunits" as well .. and maybe some future additional
similar requirements.

So if you do all that ... I have nothing against "aunits"!

regards: Janne.


Kal [b17c0de <at> gmail.com] kirjoitti: 
> Hi everyone,
> Is there any chance you will consider my +aunits= features (ticket
> #121). I am willing to donate any time required to get this integrated.
> Thanks!
> Kal
> _______________________________________________
> Proj mailing list
> Proj <at> lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj

Proj mailing list
Proj <at> lists.maptools.org

kochizufan | 18 Feb 17:52 2014

What is the difference between proj and cs2cs command?


I've believed for a long time that,

proj +proj=some_definition


cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to

is same meaning.

But I tried this:

proj -f %.6f +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001
+x_0=400000 +y_0=0 +ellps=clrk80 +towgs84=725,685,536,0,0,0,0 +units=m
-61.35 15.42 0
469737.400664	1704564.737169 0

cs2cs -f %.6f +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to
+proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0
+ellps=clrk80 +towgs84=725,685,536,0,0,0,0 +units=m +no_defs
-61.35 15.42 0
468774.391590	1704156.539473 14.452425

These cases return different results.

And, I used pj_transform C function directly, from WGS84 to target
projection, it returns same result
with proj. 
(469737.400664  1704564.737169)

What is the difference between proj and cs2cs? 
My task is making test case of pj_transform function, which result I should


View this message in context: http://osgeo-org.1560.x6.nabble.com/What-is-the-difference-between-proj-and-cs2cs-command-tp5104516.html
Sent from the PROJ.4 mailing list archive at Nabble.com.
Proj mailing list
Proj <at> lists.maptools.org

Kal | 16 Feb 21:06 2014

+aunits feature?

Hi everyone,
Is there any chance you will consider my +aunits= features (ticket
#121). I am willing to donate any time required to get this integrated.
Proj mailing list
Proj <at> lists.maptools.org

Rob Skelly | 14 Feb 02:04 2014

NAD83(Original) Transformation with Heights

Hello List,

I'm trying to convert some LiDAR data from NAD83(Original), with orthometric heights based on CGVD28, to WGS84 with ellipsoidal heights, and to NAD83(CSRS) with orthometric heights based on CGVD2013.

I've been verifying my output using NRCAN's online GPS-H tool (http://webapp.geod.nrcan.gc.ca/geod/tools-outils/gpsh.php), and not getting good results.

To begin with, I use las2las to convert my data to WGS84, using parameters:

--a_srs '+proj=utm +zone=12N +ellps=GRS80 +geoidgrids=HT2_0.gtx'
--t_srs EPSG:4326

Then I go from WGS84 to NAD83(CSRS) with the following:

--a_srs EPSG:4326
--t_srs +proj=utm +zone=12N +geoidgrid=CGG2013n83.gtx +ellps=GRS80 +datum=NAD83

Now, I don't actually have a decent way of verifying my output, other than with GPS-H, and the disagreement seems significant. For example, the original elevation at 494618.88mE, 6306815.69mN is 469.95m, and I get an ellipsoidal height at -111.088356°W, 56.90501°N of 444.52m. GPS-H claims that the orthometric height at that position (for HT2_0, the closest approximation to CGVD28) should be 470.004m, a difference of 54cm.

I realize there are a lot of points where this process could be going wrong, or that it's not going wrong at all. I wonder if anyone can give me any insight or fill in any obvious holes in my reasoning. If I've left out any pertinent information, let me know.

Note that I could just use GPS-H to process the points, but I need to automate the conversion of a lot of data, and establish a reusable process to do so. And, I'd like to use this opportunity to learn about proj!


Proj mailing list
Proj <at> lists.maptools.org
Alain Orsoni | 12 Feb 14:29 2014

Geoid models egm96 et egm08 in PROJ4


We think there could be a problem in the egm96_15.gtx and egm08_25.gtx files provided on ftp://ftp.remotesensing.org/proj/vdatum/

We noticed a difference between the values we used to get and the values provided by PROJ4 using these gtx files.
It seems that both files have a shift (although it is different for each of the files.

Probably the reason is related to the egm08togtx.py and egm96togtx.py :

=============   in egm96togtx.py the line
ds_gtx.SetGeoTransform( (-180.25,0.25,0,90.125,0,-0.25) )
should probably become
ds_gtx.SetGeoTransform( (-180.125,0.25,0,90.125,0,-0.25) )
assuming that the grid value is given for the center of a cell (0.25°x0.25° cell). 

This hypothesis seems to be confirmed by the following verification:
on the web site http://www.unavco.org/community_science/science-support/geoid/
which is a geoid height calculator, if you enter 0 0 and 0 values for long lat and height you will get a -17.162
orthometric height

If you use the following command

cs2cs +init=epsg:4326 +to +proj=longlat +datum=WGS84 +geoidgrids=egm96_15.gtx +units=m +no_defs
you will get the following result
0 0 0
0dE	0dN -17.120
while if you shift your value of -0.125
-0.125 0 0
0d7'30"W	0dN -17.160

=============  in egm08togtx.py the line
    ps_25 = 2.5 / 60.0
    ds_gtx.SetGeoTransform( (-180 - ps_25,ps_25,0,90 + ps_25,0,-1 * ps_25) )

should become

    ds_gtx.SetGeoTransform( (-180 - (ps_25/2),ps_25,0,90 + (ps_25/2),0,-1 * ps_25) )


If some one as encountered the same problem and confirm or could unconfirm this hypothesis, it would be helpful.


Proj mailing list
Proj <at> lists.maptools.org

Hamish | 9 Feb 21:59 2014

Have LCC center terms, need PROJ.4 terms


I'm working with some weather model (WRF) output and I'll like to
import the results into proj4-aware software (gdal or grass gis).

The LCC model domain specification takes center of region lat/lon as
input, not +lat_0, +lon_0. The values I've been given for +lat_1 and
+lat_2 are equal, and near to +lon_0. That doesn't look quite
right but the model seems happy enough with it:

the values I was given:

CEN_LAT = 40.26405f ;
CEN_LON = -76.99222f ;
TRULAT1 = 39.338f ;
TRUELAT2 = 39.338f ;
MOAD_CEN_LAT = 39.33799f ;
STAND_LON = -84.f ;

What I do have of value is two of the raster grids put out by the
model which contain the lat and lon of each grid point, and these
seem correct. I also have two raster outputs for convergence angle,
ie the angle between grid north and geographic north at each grid
point. I know that +x_0 and +y_0 are 0, and that +ellps is a sphere,
the WRF code lists about 4 different radii for +a, but that's a
secondary issue. Grid cell size is down to 4km, so I'm not overly
concerned with datum issues on the order of 100m.

What I'd like to do is to take the center-of-region LCC definition
and solve for the +lon_0, +lat_0, +lat_1, +lat_2 terms so I can
process the model's output grids in their native projection.

I could treat the two lat,lon grids as source data to run a manual raster
reprojection and interpolate into a proj4-know projection, but that's
computationally expensive and lossy.

So far I've gotten close by choosing the bottom-left grid point of
the lat and lon grids, and reading off those values to use for the
+lat_0 and +lon_0 terms. Then by educated-guessing I picked some
+lat_1 and +lat_2 terms and iteratively have gotten pretty close to
being able to overly a geographic grid on top of contour lines from
the lat and lon raster grids. Of course the terms I found are no
where near those I was provided with or that I can see in the model's
domain setup files.

my guess:
name       : Lambert Conformal Conic
proj       : lcc
a          : 6371229.000
es         : 0.0
f          : 0.0
lat_0      : 33.0624
lat_1      : 20
lat_2      : 45
lon_0      : -87.0188
x_0        : 0.0000
y_0        : 0.0000
unit       : meter
units      : meters
meters     : 1.0

This trial and error gets me to within about 5% (see attatched), but
I'd really like to do the job properly and calculate the exact terms.
I have access to the Fortran code used in the model to make the lat
and lon grids, could I use them to also derive the proj4 terms?

thanks for any ideas or advice on how to proceed.


ps- sorry for any crazy formatting, yahoo mail's web app is broken yet again and won't let me change to ascii or have any control over the formatting. grumble!
Proj mailing list
Proj <at> lists.maptools.org