Results of layers created with a custom CRS are loaded with a generated CRS in QGIS

Question

I have created a custom CRS for El Salvador using this parameters in a proj.4 string, based on this question:

+proj=lcc +lat_1=13.31666666666667 +lat_2=14.25 +lat_0=13.783333 +lon_0=-89 +x_0=500000 +y_0=295809.184 +k_0=0.99996704 +ellps=clrk66 +units=m +no_defs +towgs84=0,125,194,0,0,0,0 +wktext

I added +wktext to my string based on this answer. I set this CRS for my layers. I have also set the project CRS to this custom CRS. When I perform an operation using some tools, the result is loaded in the layers panel with a generated CRS. This is the generated CRS (sorry it's in spanish):

It's similar to the original, but without the +towgs84 and +wktext parameters.

It works when I save a shapefile or when I use some GDAL Tools that ask for a CRS. It doesn't work if I use some GDAL tools like GDAL Rasterize or some GRASS tools like r.surf.contour. Is there a way all the layers created with inputs with a custom CRS are loaded with this custom CRS?

Are you sure the definition is correct? EPSG:5460 matches everything except lat_1 and lat_2. Instead they says that it's a scale factor case which means lat_1 and lat_2, if required should be set to the same value as lat_0. Otherwise, use lat_1 and lat_2 but drop k_0. You're basically applying the scale factor twice. — mkennedy, Aug 28 '17 at 21:12
I downloaded a shapefile created by El Salvador National Agency. The .prj: PROJCS["IDGES_rev",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",295809.184],PARAMETER["Central_Meridian",-89.0],PARAMETER["Standard_Parallel_1",13.31666666666667],PARAMETER["Standard_Parallel_2",14.25],PARAMETER["Scale_Factor",0.99996704],PARAMETER["Latitude_Of_Origin",13.783333],UNIT["Meter",1.0]] — Ernesto561, Aug 28 '17 at 21:50
Well, I don't know. I thinl there is no a EPSG code for this particular CRS. It seems EPSG:5460 (https://epsg.io/5460-6891) is a 1SP projection. — Ernesto561, Aug 28 '17 at 21:59
LCC 1SP with a scale factor is basically the same as LCC 2SP, see http://geotiff.maptools.org/proj_list/lambert_conic_conformal_2sp.html. — AndreJ, Aug 29 '17 at 05:49
Regarding your question, the custom CRS might get lost because you are using external commands inside the QGIS GUI. You might assign the previous custom CRS to the result in a second step, unless the result gets misplaced, — AndreJ, Aug 29 '17 at 05:55
Thanks @AndreJ. It seems that in a 2SP projection the scale factor should be dropped. I'll dig a little bit more in this topic, but it seems the agency parameters should be changed. — Ernesto561, Aug 29 '17 at 11:09
@mkennedy Your are right, the EPSG:5460 matches almost everything if a 1SP is used, but the to_wgs84 parameters are different. The datum is also different (the Ocotepeque datum of 1935 ), although in this document http://www.asprs.org/a/resources/grids/07-2005-elsalvador.pdf transformations to WGS84 are mentioned, so I guess are taken as equivalents.I think your are right about the scale factor thing, and probably the agency is wrong. I have to create a custom CRS anyway.. — Ernesto561, Aug 29 '17 at 11:20
@AndreJ The workaround is changing manually the CRS in the properties as you pointed out, and that's what I've doing. So it seems there is no other way the layers are loaded with the custom CRS. Would this be a bug? — Ernesto561, Aug 29 '17 at 11:23
Dropping the datum shift is indeed a bug, which should be reported. Please include the tool you used, and some test data. — AndreJ, Aug 29 '17 at 16:07
Can you get some "official" information on Lambert NAD27? It seems this is used as a cadastral CRS in El Salvador, instead of EPSG:5460. The two differ by about 400m, due to the different datum shifts. Maybe EPSG can launch a new EPSG code for it. — AndreJ, Aug 31 '17 at 13:07
@AndreJ The "official" parameters are the ones I posted in the question. I extracted them from a prj downloaded from http://www.cnr.gob.sv/geoportal-cnr/, an official online information repository. However, it seems there is something wrong with the scale factor, as you two pointed out. I'll get in touch with them about this topic. Regarding EPSG:5460, for some reason, the towgs84 parameters are different from the ones in the reference they indicate in http://www.epsg-registry.org/. So, I think I'll submit a request to create a new EPSG code with Lambert NAD27 with 2SP. — Ernesto561, Aug 31 '17 at 15:05
For EPSG, CRS parameters and datum shift parameters are handled separately. Ocotepeque datum has a different datum shift than NAD27. So the new EPSG code should contain NAD27 as base datum. The rest is identical with EPSG:5460. Perhaps @mkennedy can speed up things integrating Lambert NAD27. — AndreJ, Aug 31 '17 at 15:11
Did you get your downloads from http://www.cnr.gob.sv/geoportal-cnr/ ? They offer WGS84 based data too. — AndreJ, Sep 01 '17 at 07:07
Thanks @AndreJ Yes, I downloaded the data from that site. i downloaded all the datasets, geographic and projected. The new EPSG code would be very useful. Anyway, I'll report a bug after double checking the issue about the custom CR in QGIS, for all the people that have similar problems. Thank you so much for your help. — Ernesto561, Sep 01 '17 at 11:45

score 1 · Accepted Answer · answered Sep 01 '17 at 16:26

Roundup answer from the comments:

It seems that datum shift and +wktext are arbitrary supported in GDAL and QGIS.

Take a screenshot in Web Mercator and warp it to your custom CRS:

gdalwarp -overwrite -dstnodata 0 -s_srs EPSG:3857 -t_srs "+proj=lcc +lat_1=13.31666666666667 +lat_2=14.25 +lat_0=13.783333 +lon_0=-89 +x_0=500000 +y_0=295809.184 +k_0=0.99996704 +ellps=clrk66 +units=m +no_defs +towgs84=0,125,194,0,0,0,0 +wktext" -of GTiff ElSalvador.png ElSalvadorNAD27.tif 
gdalsrsinfo ElSalvadorNAD27.tif >out.txt

and you get:

PROJ.4 : '+proj=lcc +lat_1=13.31666666666667 +lat_2=14.25 +lat_0=13.783333 +lon_0=-89 +x_0=500000 +y_0=295809.184 +ellps=clrk66 +towgs84=0,125,194,0,0,0,0 +units=m +no_defs '

OGC WKT :
PROJCS["unnamed",
GEOGCS["Clarke 1866",
    DATUM["unknown",
        SPHEROID["clrk66",6378206.4,294.9786982138982],
        TOWGS84[0,125,194,0,0,0,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",13.31666666666667],
PARAMETER["standard_parallel_2",14.25],
PARAMETER["latitude_of_origin",13.783333],
PARAMETER["central_meridian",-89],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",295809.184],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]

Datum shift is correctly implemented in Proj.4 string and WKT.

Take a vector file in WGS84, and convert it to the same CRS using ogr2ogr:

ogr2ogr -s_srs EPSG:4326 -t_srs "+proj=lcc +lat_1=13.31666666666667 +lat_2=14.25 +lat_0=13.783333 +lon_0=-89 +x_0=500000 +y_0=295809.184 +k_0=0.99996704 +ellps=clrk66 +units=m +no_defs +towgs84=0,125,194,0,0,0,0 +wktext" LimitesNAD27.shp dptoA_WGS_1984.shp
gdalsrsinfo LimitesNAD27.prj >>out.txt

outputs:

PROJ.4 : '+proj=lcc +lat_1=13.31666666666667 +lat_2=14.25 +lat_0=13.783333 +lon_0=-89 +x_0=500000 +y_0=295809.184 +ellps=clrk66 +units=m +no_defs '

OGC WKT :
PROJCS["Lambert_Conformal_Conic",
GEOGCS["GCS_Clarke 1866",
    DATUM["unknown",
        SPHEROID["clrk66",6378206.4,294.9786982]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",13.31666666666667],
PARAMETER["standard_parallel_2",14.25],
PARAMETER["latitude_of_origin",13.783333],
PARAMETER["central_meridian",-89],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",295809.184],
UNIT["Meter",1]]

So the datum shift gets dropped. If you do the same in QGIS, a .qpj file is created, including the datum shift values. This is good inside QGIS, but not if you run external commands.

As a sidenote, you see that the scale factor gets dropped. According to http://geotiff.maptools.org/proj_list/lambert_conic_conformal_2sp.html, LCC 1SP with a scale factor is basically the same as LCC 2SP.

If you run gdalsrsinfo on the original LambertNAD27 files from http://www.cnr.gob.sv/geoportal-cnr/, you get:

PROJ.4 : '+proj=lcc +lat_1=13.31666666666667 +lat_2=14.25 +lat_0=13.783333 +lon_0=-89 +x_0=500000 +y_0=295809.184 +datum=NAD27 +units=m +no_defs '

which does not help at all, because datum=NAD27 redirects to the grid shift files for USA and Canada, ignoring Central American states that used NAD27 as well. EPSG favours Ocotepeque datum for El Salvador, which contradicts the shapefile CRS info.

So my best advice is to use the WGS84 data from http://www.cnr.gob.sv/geoportal-cnr/ , convert them to UTM 16N (EPSG:32616), and work on with that. You will have no datum shift problem, unless you need sub-meter accuracy.

BTW the 3-parms datum shift for NAD27 is not very accurate, and Ocotepeque possibly neither. El Salvador has now introduced SIRGAS2007 for sub-meter accuracy, but no transformation parameters from old to new data have yet been published.

Thanks! It's a very detailed answer. Regarding using WGS84 with UTM, the problem is that there is a lot of information created with the Lambert 2SP CRS. My workaround right now it's just simply assign the correct CRS to the layers after every operation, but it's not practical. Besides it seems the custom CRS generates problems saving to permanent files, see this: https://issues.qgis.org/issues/17080#change-81925. Finally, if I'm understanding, ogr2ogr has a problem with preserving the CRS but GDAL doesn't, am I right? — Ernesto561, Sep 05 '17 at 12:21
Ogr2ogr is part of GDAL; it rather depends on the file format: shapefile tends to loose datum shifts, because ESRI does not save datum shifts in shapefiles. I did a follow-up question to your subject here: http://osgeo-org.1560.x6.nabble.com/gdal-dev-Which-Proj-4-transforms-are-available-in-GDAL-td5332472.html. My advice regarding UTM would be a once-and-for-all-time-reprojection, no data would get lost. — AndreJ, Sep 05 '17 at 14:31

Results of layers created with a custom CRS are loaded with a generated CRS in QGIS

1 Answers1

Linked