Crop raster from vector problem

Question

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:

Photo 2:

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

Answer 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.