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

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?

Answer

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.