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