Tuesday, 14 July 2015

gdal - Georeferencing image from Google Maps using gdal_translate?


I am trying to geo-reference an image using gdal_translate in a command prompt. Currently I have the following


gdal_translate.exe -of GTiff -a_nodata 0 --config GDAL_DATA "Config location" -a_ullr -77.057275 38.872024 -77.054672 38.869986 -a_srs EPSG:3857 "Input.tif" "Output.tif"

Where the config location is the location of the gdal data files.


The issue I am having is when I open the image it has lost its aspect ratio. I believe it has something to do with the projection but I am unsure.


gdalinfo.exe gives the following info.



Driver: GTiff/GeoTIFF
Files: Output.tif
Size is 1950, 1950
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 = (-77.057275000000004,38.872024000000003)
Pixel Size = (0.000001334871795,-0.000001045128205)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
TIFFTAG_SOFTWARE=TwelveMonkeys ImageIO TIFF writer DEV
TIFFTAG_XRESOLUTION=72
TIFFTAG_YRESOLUTION=72

Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( -77.0572750, 38.8720240) ( 0d 0' 2.49"W, 0d 0' 1.26"N)
Lower Left ( -77.0572750, 38.8699860) ( 0d 0' 2.49"W, 0d 0' 1.26"N)
Upper Right ( -77.0546720, 38.8720240) ( 0d 0' 2.49"W, 0d 0' 1.26"N)
Lower Right ( -77.0546720, 38.8699860) ( 0d 0' 2.49"W, 0d 0' 1.26"N)
Center ( -77.0559735, 38.8710050) ( 0d 0' 2.49"W, 0d 0' 1.26"N)
Band 1 Block=1950x1 Type=Byte, ColorInterp=Red
NoData Value=0

Band 2 Block=1950x1 Type=Byte, ColorInterp=Green
NoData Value=0
Band 3 Block=1950x1 Type=Byte, ColorInterp=Blue
NoData Value=0

Answer



As already pointed out in the comments, the -a_ullr parameters have to be in Google mercator projection.


So put your two coordinate pairs into a text file named wgs84.txt, and create a batch file with this content:


cs2cs +init=epsg:4326 +to +init=epsg:3857 -f "%%.2f" merc.txt

Run that batch inside the OSGeo4W shell (or Linux terminal), and you will get:



-8577976.61 4703356.65 0.00
-8577686.85 4703065.26 0.00

The third column is the height above ellipsoid, you don't need that here. Put the first two values into your gdal_translate command line, and you should get the correct georeferencing:


gdal_translate -of GTiff -a_nodata 0 -a_ullr -8577977 4703357 -8577687 4703065 -a_srs EPSG:3857 Input.tif Output.tif

If you are running the command inside the OSGeo4W shell, no config location is needed. You might have to add the full path to your source and destination file in quotes, or move to the folder with cd. Running on Linux might need different handling on that.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...