[gdal-dev] Converting from Lambert (LCC) to latlong/mercator

[Brian Murray]

I finally have positive results, that I figured I should share here
(especially since this thread is already indexed on google...)

First, I created the world file for my map. I calculated the pixel
size by getting the number of pixels between each line of longitude. I
found this to be 65.5 (1 division being 1/2 a degree, or 30 minutes,
which translates to about 55.56km. 848 pixels / 55560 = 65.5).

You may ask, how do I calculate the false easting/northing? I have no
clue. I just dumped some numbers into a website
(http://www.gazza.co.nz/datumshift.html) and viola, they were valid. I
don't even know what false easting/northings are. I gather they are
something to do with how many meters off of 0 you are, but really, no
idea.

$ cat north.wld
65.5
0
0
-65.5
8223265.173
4204567.856

Note, this world file is a companion to north.tif. While the false
east/northings are correct for any file with at the same location, the
pixels per meter measure (65) is different depending on the resolution
of the file.

Next, I wipe out all of my existing working files. gdal will _NOT_
overwrite existing files. I found this out the hard way.

rm -rf nn.* n-merc.*

Next, I converted my simple tif file to a geotif file. Photoshop added
some xres data, which I wiped, but I don't believe it is necessary.
a_ullr seems to break things, so I opted for the world file from above
instead.

echo '*** Translate ***'
gdal_translate north.tif nn.tif -of gtiff -CO TFW=YES -mo
TIFFTAG_XRESOLUTION= -mo TIFFTAG_YRESOLUTION= -a_srs "+proj=lcc
+lat_0=56 +lon_0=-115 +lat_1=49d20 +lat_2=54d40 +y_0=4204567.856
+x_0=8223265.173 +k_0=1.17230169 +datum=NAD27 +ellps=clrk80"

Note, I created a worldfile (-co TFW=YES), but this is not required.
GeoTIFF puts all of that data within the metadata of the file.

Calculating the definition (-a_srs) was rather difficult. The websites
I found were rather poorly documented. They had a lot of very
technical details, and very few laymans terms (for people like me).

lat_0 and lon_0 are the origins. That is, the lat and lon from the top
left of the image.
lat_1 and lat_2 are the standard parallels. The map I have had these
listed on it. While I cannot confirm this, these are often lines of
latitude that are 30% and 60% of the map. These are where the
"invisible cone" intersect the plane of the sphere. Everywhere between
these lines is above the cone, where all lines outside of these lines
is behind the cone. At least, thats how I understand it.
x_0 and y_0 are frustrating. I have no idea how to calculate these.
Northing is the X, and Easting is the Y, from the values from the
above website. I do not know if the website for proj is wrong, or if
the calculator is wrong. All I know is I had them backwards before,
and nothing worked right.
I don't think the k_0 is necessary, or even used, but the calculator
gave it to me, so I entered it.
Datum was written on my map as "topographic data". It was not written
with the other data, but I found it on the bottom of the map key. I
took a blind stab at the ellps.

Now I have a file called nn.tif. It contains all of the info that
properly described my map. All I had to do was tell gdalwarp what to
warp it to, and viola, out came a proper map.

echo '***   Warp    ***'
gdalwarp nn.tif n-merc.tif -t_srs "+proj=merc +datum=WGS84"

This tells gdalwarp to produce a mercator projection. Then I opened
it, and viola, it was "bent", where it was flat on the left, and
twisted all the way along to the right.

Also, you can cut out the gdal_translate, and use 1 line with the
world file instead:

gdalwarp north.tif n-merc.tif -s_srs "+proj=lcc +lat_0=56 +lon_0=-115
+lat_1=49d20 +lat_2=54d40 +y_0=4204567.856 +x_0=8223265.173
+k_0=1.17230169 +datum=NAD27 +ellps=clrk80" -t_srs "+proj=merc
+datum=WGS84"

This does the same steps, but overrides any data within the geotiff
tags. This should (and does) take a bit longer. If you have 1 map, and
want to do several projections, create an intermediary and warp
multiple times.