0

This is somewhat related to a previous question of mine.

My goal here is to clip Landsat imagery to a rectangle that I've defined by four points of latitude and longitude in WGS84. To do this, I've used this python module to convert the points of the rectangle into utm coordinates suitable for use as arguments to gdalwarp. My question is how do I determine what the correct EPSG code is for gdalwarp? I assumed it would be a simply mapping of UTM Zone to EPSG (for example, coordinates in UTM Zone 11N (southern california) would map to http://www.spatialreference.org/ref/epsg/wgs-84-utm-zone-11n/). However, when I try to use this code in my gdalwarp command, it hangs trying to create an image much too large. When I change the EPSG code to http://spatialreference.org/ref/epsg/nad83-utm-zone-11n/ it works and gives me the clipped imagery I expect. What I don't understand is that the info page on the USGS web site implies that 32611 is the correct EPSG code for this specific landsat tile, but it seems that 26911 is the actual correct EPSG code. And further, why don't either of these EPSG codes appear in the output from gdalinfo for the file?

The command I'm using to clip the images is:

gdalwarp -of gtiff -t_srs EPSG:<EPSG Code goes here> -te <bounds in utm>  input.tif output.tif

The output from gdalinfo for my input image is:

Size is 15381, 15681
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 = (-13185978.181983433663845,4059663.631883424706757)
Pixel Size = (18.129383469292915,-18.023985075689048)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-13185978.182, 4059663.632) (118d27' 5.97"W, 34d13'46.30"N)
Lower Left  (-13185978.182, 3777029.522) (118d27' 5.97"W, 32d 6'16.03"N)
Upper Right (-12907130.135, 4059663.632) (115d56'48.20"W, 34d13'46.30"N)
Lower Right (-12907130.135, 3777029.522) (115d56'48.20"W, 32d 6'16.03"N)
Center      (-13046554.158, 3918346.577) (117d11'57.08"W, 33d10'24.35"N)
Band 1 Block=15381x1 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=15381x1 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=15381x1 Type=Byte, ColorInterp=Blue
  NoData Value=0

Sorry if this is a dumb or obvious question, I'm still new to GIS.

[EDIT]

Finally figured it out!

I was using landsat-utils to process the landsat imagery down to a single rgb tiff. In the process, the image was getting mapped to web mercator from the UTM of the original input imagery. I'm still not sure why I was getting reasonable results with UTM, but by using this gdalwarp command, everything comes out more or less as expected.

gdalwarp -of gtiff -t_srs EPSG:3857 -te <bounds in web mercator>  input.tif output.tif
  • Looks like landsat-util only supports EPSG:3857 as the output CRS for processed imagery. Which is useless other than for web mapping, but at least alternative CRSs are on their TODO list... As an alternative workflow, you can use the --clip argument and pass lon/lat coordinates in WGS84 – user2856 Mar 1 '16 at 22:10
1

Your gdalinfo output shows that the input image is in EPSG:3857 aka Web Mercator which is mostly used by online mapping apps, i.e Google Maps.

I don't know why your gdalwarp command works as expected with EPSG:26911, something else may be going on. It might help if you edited your question to specify where you got the Landsat image from and what (if any) processing you have done to it as the number of bands and the byte data type show that it's not raw Landsat data.

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.