4

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):

enter image description here

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?

  • 1
    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
  • Thank you! So, is EPSG confused or the agency?? Hmmm. – mkennedy Aug 28 '17 at 21:53
  • 1
    LCC 1SP with a scale factor is basically the same as LCC 2SP, see geotiff.maptools.org/proj_list/lambert_conic_conformal_2sp.html. – AndreJ Aug 29 '17 at 5:49
  • 1
    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
1

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: 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 5 '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: osgeo-org.1560.x6.nabble.com/…. My advice regarding UTM would be a once-and-for-all-time-reprojection, no data would get lost. – AndreJ Sep 5 '17 at 14:31

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.