1

I've been having a problem with GDAL/OGR lately mixing up projections EPSG:3857 and EPSG:54004. I'm cropping some 3857-projected raster files (reprojected using gdalwarp from UTM) by 3857-projected vector files (reprojected using ogr2ogr from KML files - EGPS:4326). I checked gdalinfo and ogrinfo for both and they seem to be the right projection. The result is a clipped raster that only contains part of what I want due to the offset between 3857 (which it seems the raster stays in) and 54004 (which it seems to think the vector is in). So generally only a part of (some of the bottom) of the desired raster will show due to the 54004 vector being offset north of where the 3857 raster is.

I did see a similar problem in this post, but this made me think it was a QGIS issue as opposed to underlying GDAL/OGR which I've exclusively used for my data processing. While the workaround suggested there via QGIS is effective, I will be doing this with many files, so it seems impractical for my requirements. When I check GDAL version I'm on: "GDAL 1.11.3, released 2015/09/16". I installed this using homebrew (brew install gdal) on a mac, and this seems to be the latest version available via this method from kyngchaos.

So: Is it possible I'm doing something wrong with my projecting/processing (see info below) or is it perhaps my version of GDAL which may be the problem as seen here? If it is the gdal version, please advise the best way to upgrade this on a mac (since home-brew insists I'm using the latest version with 'brew upgrade gdal'?

Here is my clipping command:

gdalwarp -crop_to_cutline -cutline vector.shp raster.vrt output.tif

Also below are example 'gdalinfo' and 'ogrinfo -so' outputs, respectively.

RASTER:
Driver: VRT/Virtual Raster
Files: raster.vrt
Size is 63319, 54925
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-11716505.067412510514259,4866089.724105604924262)
Pixel Size = (1.323481832251734,-1.323481832251734)
Corner Coordinates:
Upper Left  (-11716505.067, 4866089.724) (105d15' 4.16"W, 40d 0' 3.65"N)
Lower Left  (-11716505.067, 4793397.484) (105d15' 4.16"W, 39d29'56.26"N)
Upper Right (-11632703.521, 4866089.724) (104d29'54.07"W, 40d 0' 3.65"N)
Lower Right (-11632703.521, 4793397.484) (104d29'54.07"W, 39d29'56.26"N)
Center      (-11674604.294, 4829743.604) (104d52'29.12"W, 39d45' 1.60"N)
Band 1 Block=128x128 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=128x128 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=128x128 Type=Byte, ColorInterp=Blue
 NoData Value=0



VECTOR:
INFO: Open of `vector.shp'
  using driver `ESRI Shapefile' successful.

Layer name: vector
Geometry: Polygon
Feature Count: 1
Extent: (-11705558.991470, 4803872.365048) - (-11639276.914951, 4859716.734835)
Layer SRS WKT:
PROJCS["WGS_84_Pseudo_Mercator",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1],
    PARAMETER["latitude_of_origin",0.0]]
Name: String (80.0)
Descriptio: String (80.0)

Output from 'gdalsrsinfo -e vector.shp':
Warning 1: EPSG detection is experimental and requires new data files (see bug #4345)

EPSG:-1

PROJ.4 : '+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs '

OGC WKT :
PROJCS["WGS_84_Pseudo_Mercator",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1],
    PARAMETER["latitude_of_origin",0.0]]
Improve this question
9
  • 1
    I'm having a problem editing the question because the ProjCRS definition is neither Esri:54004 (note that it's not EPSG) nor EPSG:3395 which it's closer to except for the name. Notice that the raster version has an EXTENSION definition which specifies a sphere. The overall ProjCRS definition doesn't which is why you're getting mismatched results.
    – mkennedy
    Commented Sep 1, 2016 at 22:17
  • Thanks for your comment, but I'm not sure where to go with it - I'm not sure what the EXTENSION on the raster means, but it seems the raster behaves correctly (e.g. opens in QGIS as a EPSG:3857.) The culprit would appear to be the vector shape file which opens in QGIS as EPSG:54004 (despite the info shown in the vector ogrinfo above). I've used gdalwarp and ogr2ogr both with the '-t_srs EPSG:3857' argument to obtain the raster and vector, respectively. Have I done anything wrong by doing this?
    – GISBeginner
    Commented Sep 1, 2016 at 22:45
  • For the shapefile, either edit the .prj file and insert the EXTENSION bit or copy the whole version from the raster and replace the .prj contents (keep it was one line). It's perfectly valid to have a Mercator ProjCRS on an ellipsoid which is what the vector's .prj is saying that it's using.
    – mkennedy
    Commented Sep 1, 2016 at 22:49
  • Tried with no luck-but anyway, I'm not sure we're on the right track here: when I open the shapefile (without editing the .prj file), it appears as EPSG:54004 in QGIS. If I then view shapefile properties in QGIS and change the CRS to 3857, click okay and then save this to a new file, I get the exact same ogrinfo as shown above (there is, however, a .qpj file now associated with the shapefile). Using this new shapefile, however, the clipping works fine. So this workaround works (as in linked post above), but I don't intend to do this manually for the many many shape files I have to process.
    – GISBeginner
    Commented Sep 1, 2016 at 23:13
  • You could copy the .qpj file (since that seems to work) to the other shapefiles. I don't know QGIS that well so can't make a better suggestion.
    – mkennedy
    Commented Sep 1, 2016 at 23:17

1 Answer 1

Reset to default
1

You wrote that you are using GDAL 1.11.3. GDAL version 2.2 (master) gives different results. I made a test shapefile with ogr2ogr by using parameter -a_srs epsg:3857. Now the gdalsrsinfo report looks like this:

gdalsrsinfo -e test.shp
Warning 1: EPSG detection is experimental and requires new data files (see bug #
4345)

EPSG:-1

PROJ.4 : '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=
0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs'

OGC WKT :
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +
x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]

I believe that your problems will disappear once you can update your command line GDAL into 2.x version and your QGIS is compiled with GDAL 2.x as well.

Improve this answer
1
  • Thanks user30184 - this does appear to confirm my problem, so I've marked it as the answer. I've successfully downloaded and installed GDAL 2.1.0, and I now see something similar with ogrinfo as what you're pasted above. So it seems it was a bug with OGR. However, I'm having a hard time figuring out to install QGIS against this new version of GDAL. Is this the correct process in this link: htmlpreview.github.io/?https://raw.github.com/qgis/QGIS/master/… ? Any other guidance for how to compile QGIS with GDAL 2.x for Mac OS X? Thanks for your help.
    – GISBeginner
    Commented Sep 2, 2016 at 17:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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