4

I am trying to clip a county DEM file using the PLSS Section shape files as mask layers... that is, I want it clipped to sections so that I can load my DEM into QGIS one section at a time, without having to load the whole 4GB raster for any little project I want to work on.

I have tried to do this by going to Raster > Extraction > Clipper ... however it returns an error:

ERROR 6: No translation for Lambert_Conformal_Conic to PROJ.4 format is known. ERROR 6: No translation for Lambert_Conformal_Conic to PROJ.4 format is known.

I have found a way to rasterize the shape files and use the raster calculator, however this is more of a pain and I want to understand why this is not working for me.

I assume it has something to do with the fact that I'm using a custom CRS for my county? CRS databases don't recognize it.

+proj=lcc +lat_1=43.57503293972223 +lat_0=43.57503293972223 +lon_0=-90.78333333333333 +k_0=1.0000408158 +x_0=222504.44500889 +y_0=47532.0603505207 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs

Running QGIS 2.0.1 on Macbook Air using Mavericks 10.9


I ran gdalsrsinfo on the raster file and the command line returned the same error, plus the WKT:

ERROR 6: No translation for Lambert_Conformal_Conic to PROJ.4 format is known.
PROJ.4 : ''

OGC WKT :
PROJCS["NAD_1983_HARN_WISCRS_Vernon_County_Feet",
    GEOGCS["GCS_North_American_1983_HARN",
        DATUM["D_North_American_1983_HARN",
            SPHEROID["GRS_1980",6378137.0,298.257222101]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Lambert_Conformal_Conic"],
    PARAMETER["False_Easting",730000.0],
    PARAMETER["False_Northing",155944.768],
    PARAMETER["Central_Meridian",-90.78333333333333],
    PARAMETER["Standard_Parallel_1",43.57503293972223],
    PARAMETER["Scale_Factor",1.0000408158],
    PARAMETER["Latitude_Of_Origin",43.57503293972223],
    UNIT["Foot_US",0.30480060960121924],
    VERTCS["NAVD_1988",
        VDATUM["North_American_Vertical_Datum_1988"],
        PARAMETER["Vertical_Shift",0.0],
        PARAMETER["Direction",1.0],
        UNIT["Foot_US",0.3048006096012192]]]
3

This seems to be a known bug in GDAL working in the background of QGIS about different ways of definition for LCC 1-SP and 2-SP projections:

See

Adding ESRI:: in a gdalwarp command line as described should solve the problem. Or change the shapefiles .prj to

PROJCS["NAD_1983_HARN_WISCRS_Vernon_County_Feet",
GEOGCS["GCS_North_American_1983_HARN",
    DATUM["D_North_American_1983_HARN",
            SPHEROID["GRS80",6378137,298.257222101],
            TOWGS84[0,0,0,0,0,0,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_1SP"],
    PARAMETER["latitude_of_origin",43.57503293972223],
    PARAMETER["central_meridian",-90.78333333333333],
    PARAMETER["scale_factor",1.0000408158],
    PARAMETER["false_easting",730000],
    PARAMETER["false_northing",155944.768],
    UNIT["Foot_US",0.3048006096012192]]

and run gdal_translate -a_srs new.prj src_dataset dst_dataset on the raster file.

Note that false Easting and Northing are in feet for the .prj file, and always in meters for the +proj string.


UPDATE

I found the DEM at http://relief.ersc.wisc.edu/wisconsinview/form.php and PLSS data at http://www.geocommunicator.gov/GeoComm/lsis_home/home/PLSS_download_WIcounty.htm (in NAD 83)

The PLSS is also imprinted on USGS Topo maps. So I took one township (T13N R02W), and reprojected the shapefile to Vernon WCCS with the parameters given above. The DEM comes in adf format, so I translated that to Vernon WCCS as well:

gdal_translate -a_srs "+proj=lcc +lat_1=43.57503293972223 +lat_0=43.57503293972223 +lon_0=-90.78333333333333 +k_0=1.0000408158 +x_0=222504.44500889 +y_0=47532.0603505207 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs" -of GTiff D:/Karten/shp/Laender/USA/vernon/w001000.adf D:/Karten/shp/Laender/USA/Vernon-DEM.tif

gdalwarp -q -cutline D:/Karten/shp/Laender/USA/T13NR02W.shp -crop_to_cutline -of GTiff D:/Karten/shp/Laender/USA/Vernon-DEM.tif D:/Karten/shp/Laender/USA/Vernon-DEM-clipped.tif

These are actually the command lines generated by QGIS. The result fits to the township extent:

enter image description here

Background is the La Farge Topo from 1965.

  • Thanks for your help AndreJ. I am still learning the ropes with QGIS, I have only dabbled in the command line side of things. I updated the shapefile's .prj as you suggested, but the problem did not resolve. Then I tried running gdalwarp in terminal, with ESRI:: prepended to the filename as suggested, and it is not returning errors but it does not seem to be producing any output. I tried: gdalwarp -q -cutline ESRI::"/Users/shapefile.shp" -crop_to_cutline -of GTiff "/Users/demraster.tif" /Users/output – corvus Jan 6 '15 at 15:52
  • For the shapefile, you have to delete the .qpj file QGIS creates. Otherwise it will not look into the exchanged .prj file. You have to run gdal_translate on the raster file to change from lcc to lcc_1SP (which seems to be the problem). Can you give download links to the data? I have not found the shapefiles you mention. I guess the DEM is from sco.wisc.edu/news/… – AndreJ Jan 6 '15 at 16:27
  • Thanks again, I will give that a try. DEM is a 1M LiDAR raster from WisconsinView… Vernon County is the one. I got the PLSS Sections from a buddy of mine who is a GIS professional. Not sure of the sourcing, it is a simple polygon shape file for our county. I can try to upload it somewhere if you want to have a look at it – corvus Jan 6 '15 at 17:15
  • Am I right in thinking this would be accomplished by gdal_translate -a_srs "SRS" source.tif output.tif?? what should I substitute for "SRS" to get LCC 1SP as my output? – corvus Jan 6 '15 at 17:46
  • Wow, thanks so much! Worked like a charm. I have been trying to figure this one out for so long, you have no idea how much I appreciate all of your effort! – corvus Jan 7 '15 at 4: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.