Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a png source image that's in Lambert Azimuthal Equal Area and I have the long/lat extent. So I ran gdal_translate to create tie points and convert to gtif. Then I ran gdalwarp to reproject to Web Mercator. gdalinfo shows that the bounding coordinates are ok for the gdal_translate output, but not for the gdalwarp output. Que pasa?

Here's my gdal_translate operation on the source file:

gdal_translate -of GTiff -a_srs "+proj=laea +lat_0=40. +lon_0=-105 +x_0=0 +y_0=0 +a=6370997 +rf=298.257" ^
-gcp 0 0 -119.344674 47.148697 ^
-gcp 525 0 -90.655326 47.148697 ^
-gcp 525 431 -93.609897 31.438244 ^
-gcp 0 431 -116.390103 31.438244 ^
-a_ullr -119.344674 47.148697 -93.609897 31.438244 ^
src.png tmp.gtif

(I included the bounding box long/lats for good measure but it seems to make no difference in the outcome of gdalwarp.)

Here's my gdalwarp operation on the gdal_translate output (tmp.gtif):

gdalwarp -of GTiff -s_srs "+proj=laea +lat_0=40. +lon_0=-105 +x_0=0 +y_0=0 +a=6370997 +rf=298.257" -t_srs EPSG:3857 tmp.gtif final.gtif

Here's the metadata from 'gdalinfo tmp.gtif':

Driver: GTiff/GeoTIFF
Files: tmp.gtif
Size is 525, 431
Coordinate System is `'
Origin = (-119.344674000000000,47.148696999999999)
Pixel Size = (0.049018622857143,-0.036451167053364)
Image Structure Metadata:
  INTERLEAVE=PIXEL
**Corner Coordinates:
Upper Left  (-119.3446740,  47.1486970) 
Lower Left  (-119.3446740,  31.4382440) 
Upper Right ( -93.6098970,  47.1486970) 
Lower Right ( -93.6098970,  31.4382440) 
Center      (-106.4772855,  39.2934705)** 
Band 1 Block=525x3 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA 
Band 2 Block=525x3 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA 
Band 3 Block=525x3 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA 
Band 4 Block=525x3 Type=Byte, ColorInterp=Alpha

Corner coordinates look good. But here's the metadata from 'gdalinfo final.gtif':

Driver: GTiff/GeoTIFF
Files: final.gtif
Size is 579, 355
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 = (-11688702.286414048000000,4866004.054033085700000)
Pixel Size = (0.057992644981843,-0.057992644981843)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
**Corner Coordinates:
Upper Left  (-11688702.286, 4866004.054) (105d 0' 5.04"W, 40d 0' 1.53"N)
Lower Left  (-11688702.286, 4865983.467) (105d 0' 5.04"W, 40d 0' 1.02"N)
Upper Right (-11688668.709, 4866004.054) (105d 0' 3.95"W, 40d 0' 1.53"N)
Lower Right (-11688668.709, 4865983.467) (105d 0' 3.95"W, 40d 0' 1.02"N)
Center      (-11688685.498, 4865993.760) (105d 0' 4.49"W, 40d 0' 1.28"N)**
Band 1 Block=579x3 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA 
Band 2 Block=579x3 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA 
Band 3 Block=579x3 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA 
Band 4 Block=579x3 Type=Byte, ColorInterp=Alpha

Also notice that the pixel size is suspect in the final image - both dx and dy are the same.

share|improve this question
    
I'm suspicious about the use of +a and +rf. 6370997 is a sphere (almost authalic based on Clarke 1866), not an ellipsoid. So you might try +b=6370997, not +rf=298.257. A lot of US government data in Lambert azimuthal uses a sphere-only implementation. –  mkennedy Aug 5 '11 at 23:31
    
Thanks for the feedback, Mike. I'll try that and compare. –  Jeffo Aug 6 '11 at 22:35
add comment

1 Answer 1

up vote 1 down vote accepted

The reason you seem to be getting very small values is because your original image is projected, so it uses linear units - which in this case PROJ4 I think defaults to metres - rather than degrees. So you're specifying a ~680 square metre part of the western US.

Now, how you can determine the right coordinates to use is a tricky matter because your image covers such a large area, you'll get a lot of distortion at the edges.

Sorry this isn't much of a solution - anyone else care to take over? If I have a flash of inspiration, I'll post it here.

share|improve this answer
    
Ah, that makes sense - that Proj would assume meters since it's projected. –  Jeffo Aug 6 '11 at 22:35
    
I started with a projected file and was given a center lat/lon and the spatial extent in meters. So I used proj -I (inverse) to determine lat/lon values of the corners from the initial info, and I used those lat/lon values for gcps. Could I then use proj to get northing and easting for the corners from the lat/lons and use those linear values for gcps? Or is that going to produce the distortion you mention? Thanks for looking into this! –  Jeffo Aug 6 '11 at 22:45
    
Ah, but that's circular - it will give me linear values with respect to the center of the image, which I already have. You're saying those linear values will be distorted? –  Jeffo Aug 6 '11 at 22:56
1  
I just tried it and it worked - I used the corner values in meters from the spatial extent of the original file. Indeed there is some distortion near the edges, but I think it's acceptable for this application. Sincere thanks for the answer. –  Jeffo Aug 6 '11 at 23:33
    
Ah, I'm glad you've got it sorted. stackexchange is all about giving the questioner the tools to find their own solutions... well, that's what I'm telling myself :) If you could mark the answer as accepted, that'd be cool. –  MerseyViking Aug 7 '11 at 12:39
show 1 more comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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