5

My problem:

I'm currently trying to use gdal to write GeoTIFF files. I'm however having problems with using a proj4 string to describe the projection. ("+proj=lcc +k_0=0.001 +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.333333 +ellps=bessel +x_0=400 +y_0=400 +no_defs"). Essentially the scaling factor "+k_0" is ignored by gdal. I know that instead of +k_0, I can use dataset->SetGeoTransform(adfGeoTransform) (in C++) and set this value manually. However I am surprised that doing this via the proj4 string does not seem to be possible.

If you try the following in python:

Example:

from gdal import osr
p = osr.SpatialReference()
p.ImportFromProj4("+proj=lcc +k_0=0.001 +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.333333 +ellps=bessel +x_0=400 +y_0=400 +no_defs")
print p.ExportToProj4()

the output simply drops the +k_0 parameter:

+proj=lcc +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.333333 +x_0=400 +y_0=400 +ellps=bessel +units=m +no_defs

This results in completely different projected coordinates however and makes it difficult to rely onto.

Improve this question

2 Answers 2

Reset to default
6

The value for +k_0 should be around 1.0. See the table at Snyder, Map projections, page 110. You may find the source code for lcc at https://trac.osgeo.org/proj/browser/trunk/proj/src/PJ_lcc.c

If you intend to use units different from meters, use +units or +to_meter instead.

Testing with cs2cs on 1000 1000:

echo k=1 m >out.txt
cs2cs +proj=lcc +k_0=1 +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.333333 +ellps=bessel +x_0=400 +y_0=400 +no_defs +to +init=epsg:4326 -f "%%.9f" lcc-in.txt >>out.txt
echo k=0.001 m >>out.txt
cs2cs +proj=lcc +k_0=0.001 +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.333333 +ellps=bessel +x_0=400 +y_0=400 +no_defs +to +init=epsg:4326 -f "%%.9f" lcc-in.txt >>out.txt
echo k=1 mm >>out.txt
cs2cs +proj=lcc +k_0=1 +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.333333 +ellps=bessel +x_0=400 +y_0=400 +units=mm +no_defs +to +init=epsg:4326 -f "%%.9f" lcc-in.txt >>out.txt
echo k=1 to-meter >>out.txt
cs2cs +proj=lcc +k_0=1 +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.333333 +ellps=bessel +x_0=400 +y_0=400 +to_meter=0.001 +no_defs +to +init=epsg:4326 -f "%%.9f" lcc-in.txt >>out.txt
echo k=0.9 m >>out.txt
cs2cs +proj=lcc +k_0=0.9 +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.333333 +ellps=bessel +x_0=400 +y_0=400 +no_defs +to +init=epsg:4326 -f "%%.9f" lcc-in.txt >>out.txt

delivering:

k=1 m 
13.341301026    47.505398796 0.000000000
k=0.001 m 
22.172450742    52.582612826 0.000000000
k=1 mm 
13.328035167    47.496409491 0.000000000
k=1 to-meter 
13.328035167    47.496409491 0.000000000
k=0.9 m 
13.342186463    47.505998628 0.000000000

The EPSG Guidance 7-2 lists a 2SP Michigan variant of lcc, using a scale factor and 2 parallels. The factor seems to scale the ellipsoid. In the proj definitions for EPSG 4268 and 26811 to 26813, this is been done by scaling a and b (ignoring the NAD27 datum shift to WGS84). Using +k_0=1.0000382 on the original ellipsoid returns the same values, as given in the Guidance.

See http://lists.maptools.org/pipermail/proj/2015-April/007108.html for the discussion on that.

Improve this answer
6
  • Ok. It seems the answer to my question is really: "Not every string that works in proj.4 is valid in gdal." The +to_meter seems to be what I was actually looking for when trying out +k
    – knedlsepp
    Sep 19, 2016 at 10:18
  • Why do you think that it works in proj.4? I assume wrong parameters are dropped silently.
    – AndreJ
    Sep 19, 2016 at 10:19
  • Because if I transform some coordinates (with the proj.4 api) using the string with "+k=0.001", I get the same coordinates as without, only scaled by a factor of 1e3. So adding +k=0.001 changes the result, which is what I was expecting to happen also in gdal, but there both with and without yield the same coordinates.
    – knedlsepp
    Sep 19, 2016 at 10:23
  • If I use Proj.4 cs2cs to convert from your lcc to EPSG:4326, k=0.001 and +to_meter=0.001 deliver different coordinates. The first one is totally off. See my extended answer.
    – AndreJ
    Sep 19, 2016 at 12:19
  • I was actually comparing the following two definitions: with +k: "+proj=lcc +lat_1=46N +lat_2=49N +lat_0=47.5N +lon_0=13.333333 +ellps=bessel +x_0=400 +y_0=400 +k=0.001 +no_defs" and with +to_meters "+proj=lcc +lat_1=46N +lat_2=49N +lat_0=47.5N +lon_0=13.333333 +ellps=bessel +x_0=400000 +y_0=400000 +to_meters=1000 +no_defs". Using proj4 they are equivalent. Using gdal they are not (because +k is not supported.) I agree with you that units is the better way to solve this problem.
    – knedlsepp
    Sep 19, 2016 at 13:01
3

By reading the documentation https://trac.osgeo.org/gdal/wiki/NetCDF_ProjectionTestingStatus and looking at the note about lcc it seems that +k_0 is not supported and defaults to 1 always.

Lambert conformal 1SP / 2SP

The 1SP variant corresponds to EPSG code 9801 - Lambert Conic Conformal (1SP), with CF latitude_of_projection_origin=standard_parallel and WKT scale_factor=1. A scale factor less than 1 means that there are 2 standard parallels (see www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html), but it cannot be translated to the CF 1SP variant, therefore the 2SP variant should be used instead. The 2SP variant corresponds to EPSG code 9802 - Lambert Conic Conformal (2SP).

BTW. Links to remotesensing.org do not work, someone else has taken the domain name. Use web.archive.org instead.

https://web.archive.org/web/20160305172849/http://remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html

https://web.archive.org/web/20160305140706/http://remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp.html.

Improve this answer
7
  • 1
    So this means that not every valid proj.4 string can be used as a gdal-OGRSpatialReference?
    – knedlsepp
    Sep 19, 2016 at 8:49
  • Ok. I was expecting +k to work, because if you use it in proj.4 alone, it does correctly scale using the given value. However it seems that even though proj.4 does not a have a problem combining LCC-2SP and +k, gdal does.
    – knedlsepp
    Sep 19, 2016 at 10:15
  • Probably you understand it better, but by reading the links above I concluded that 2SP does not support +k_0 at all, and if +k_0 is used with 1SP it means effectively a conversion into 2SP with two properly set standard parallels.
    – user30184
    Sep 19, 2016 at 10:36
  • Yes, the links describe gdal's behaviour when using proj4 strings. You are right that what I was trying to do is simply not supported by gdal. What threw me off however is, that if you use the proj.4 API directly, the parameter +k is supported.
    – knedlsepp
    Sep 19, 2016 at 10:48
  • But even if it is supported, is is the correct way for converting millimeters to meters if that was your real need?
    – user30184
    Sep 19, 2016 at 11:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.