3

So I have a shapefile of a vector for Philadelphia's city lines and a raster image of Philadelphia county that I got from the lovely people from USDA:NRCS. I converted the raster from a sid file to TIFF (using gdal_translate) and then reprojected the newly converted raster TIFF to 3857 (using gdalwarp).

This new raster now has the same native projection as the vector for the city limits: EPSG:3857. So I load both the shapefile for the city limits and this new raster tif that I prepped into QGIS and I fire up that useful clipper tool (Raster->Extraction->Clipper) and I go about my business of setting the raster as the source and selecting the vector as my mask layer and press "OK". #thatwassimple. The computer starts chugging along doing its thing and I'm getting pumped about my new raster coming in.

Ding! computer finished... only South Philly and a tiny part of Central Philly are extracted out--West, North, and East Philly seem to have ducked out of the party and it's all black.

Hm... that's weird, maybe I can use GDAL command line to force this, because you know, if it's done on the command line, I am a real hacker and it's not like QGIS is already doing this in the background, but I digress.

So I cmd those mofos (using gdalwarp again!) 0...10....100! DING! Alright, alright, alright, as I McConaughey myself over to the directory that I outputted the clipped TIFF. Open that sucker up and I'm with the exact same output of Southside and Central only left to represent.

I've attached two images, the first shows a rose colored polygon that I want to cut out of the raster, and the second shows the output both times.

Photo 1: rose colored polygon that I want to cut out of the raster

Photo 2: output both attempts to clip

SOURCE FILES

Shapefile used (in zip)

TIFF used (in zip)

UPDATE: I noticed that the clipped areas are actually in fact regions outside of the city and just cuts north of the city. (Note south philly no longer has Schuylkill River going through it). I tried making a polygon of this so I could upload an image for clarity, but I couldn't figure it out. Please just replicate the raster clipping with the files provided above and you too can see what I'm talking about.

This leads me to believe that there must be a projection problem, however, after I gdalinfo-ed the raster and consulting the .prj of the city_limits file, I've confirmed that they have the same aforementioned projection. #sayhwatt

  • what was the commandline statement that you used? did you use the -cutline flag? – user1269942 Dec 29 '15 at 23:49
  • Yes @user1269942, I followed the following syntax: inline gdalwarp -cutline city_limits.shp -crop_to_cutline philly_3857.tif philly_cropped.tif – Jonathan Dec 30 '15 at 2:00
  • try checking topology of the vector layer – Maxim Dubinin Dec 30 '15 at 4:25
  • Can you upload the shapefile somehwere? – AndreJ Dec 30 '15 at 6:46
  • @AndreJ here is the link: filedropper.com/citylimits – Jonathan Dec 30 '15 at 12:05
1

There seems to be a small but significant difference between the projections that makes gdalwarp stumble.

Running gdalsrsinfo philly_3857.tif > city_limits.prj and removing the first lines from it did the trick for me.

For comparison:

philly_3857.tif:

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"]]

city_limits.prj (the original):

PROJCS["WGS_84_Pseudo_Mercator",
GEOGCS["GCS_WGS_1984",
DATUM["D_WGS_1984",
SPHEROID["WGS_1984",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["Degree",0.017453292519943295]],
PROJECTION["Mercator"],
PARAMETER["central_meridian",0],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["Meter",1],
PARAMETER["standard_parallel_1",0.0]]

For a normal user, they seem to be the same. But if you delete the .q** files from the zip, the shapefile gets loaded into QGIS as EPSG:54004, and it will not align with the raster. gdalwarp does not look into the .q** files, and so it produces garbage.

  • Ok just to clarify, although defined to be the same projection, because city_limit used the nomenclature of "WGS_84_Pseudo_Mercator" whereas philly_3857 used "WGS 84 / Pseudo-Mercator", this prevented the proper clipping of the desired raster. And although it looked fine and dandy on QGIS, it failed because gdalwarp didnt look into QGIS's proprietary projection file (qpj) but was rather looking for ESRI's projection file (which file extension)? – Jonathan Dec 30 '15 at 19:20
  • I have not checked it further, but the EXTENSION parameter should make it clear to gdalwarp that the special google mercator hack of ellipsoid and sphere should be used. The .qpj contains the same, while ESRI and gdalwarp does not use it. – AndreJ Dec 31 '15 at 6:13

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.