Helper file to create a geotiff from a given image or from an automatically downloaded satellite image from google maps which is centered around a given coordinate .

create_geotiff.py

"""
BST License 
# Description
File to create geotiffs depends on gdal installed on system. Can be used under two modes:
    1) Download a satellite overview screenshot from google maps around a center
    coordinate. Can give zoom (max is usually 20, but depends on area - will return
    an image with error message if zoom is too high) and image size
    (defaults to max 640x640)
    2) Convert an existing image into a geotiff by giving UL and LR corner
    coordinates of image.
Example:
    Create geotiff from google around Barbados:
        python create_geotiff.py --CLL "13.159717, -59.657345" -z 12
    Create geotiff from existing map image around Barbados:
        python create_geotiff.py -s ~/Desktop/map_image.png --ULLR "13.1601374805, -59.6577740797, 13.1592965187, -59.6569159218"
# Author: Johanna Hansen <[email protected]>
# Edits:
    created 19 Dec 2016
# Install (tested on Python 2.7)
    pip install scipy LatLon urllib
"""
import os
from scipy.misc import imread, imsave
import math

def download_screenshot(img_path, LC, s_size_w=8000, s_size_h=8000, zoom=20):
    import urllib
    from io import BytesIO
    import LatLon
    latcenter = LC[0].strip()
    loncenter = LC[1].strip()
    # max zoom is 21, 1 is the whole earth
    url = "http://maps.googleapis.com/maps/api/staticmap?maptype=satellite&center=%s,%s&size=%sx%s&zoom=%s" %(latcenter,loncenter,s_size_w, s_size_h, zoom)
    buffer = BytesIO(urllib.urlopen(url).read())
    img = imread(buffer)
    # determine how many meters are represented in each pixel
    km_per_pixel = (156543.03392 * math.cos((float(latcenter) * math.pi) / 180.0)) / (2.0**int(zoom)) *1e-3
    h, w, _ = img.shape
    hm = h/2
    wm = w/2
    imsave(img_path, img)
    off_center_wkm = w/2.0 * km_per_pixel
    off_center_hkm = h/2.0 * km_per_pixel
    CLat = LatLon.Latitude(latcenter)
    CLon = LatLon.Longitude(loncenter)
    center = LatLon.LatLon(CLat, CLon)
    # North is 0
    # West is 270
    CT = center.offset(0, off_center_hkm)
    CB = center.offset(180, off_center_hkm)
    UL = CT.offset(270, off_center_wkm).to_string('D')
    LR = CB.offset(90, off_center_wkm).to_string('D')
    print('Calculated UL, LR coordinates of image is:', UL, LR)
    return UL, LR

def create_geotiff(img_path, UL, LR):
    if not os.path.exists(img_path):
        print("ERROR: image path to convert to geotiff does not exists: %s" %img_path)

    geotiff_path = ''.join(img_path.split('.')[:-1]) + '_geo.tiff'
    print("Creating geotiff named: %s" %geotiff_path)
    # -a_ullr ulx uly lrx lry:
    #Assign/override the georeferenced bounds of the output file.
    #This assigns georeferenced bounds to the output file, ignoring what would have been derived from the source file.
    #So this does not cause reprojection to the specified SRS.
    cmd = "gdal_translate -a_nodata 0 -of GTiff -a_srs EPSG:4326 -a_ullr %s %s %s %s  %s %s" %(UL[1], UL[0],
                                                                 LR[1], LR[0],
                                                                 img_path, geotiff_path)

    print("Running gdal command:", cmd)
    os.system(cmd)

if __name__ == '__main__':
    import argparse

    parser = argparse.ArgumentParser(
        description="""
File to create geotiffs. Can be used under two modes:
    1) Download a satellite overview screenshot from google maps around a center
    coordinate. Can give zoom (max is usually 20, but depends on area - will return
    an image with error message if zoom is too high) and image size
    (defaults to max 640x640)
    2) Convert an existing image into a geotiff by giving UL and LR corner
    coordinates of image.
Example:
    Create geotiff from google around Barbados:
        python create_geotiff.py --CLL "13.159717, -59.657345" -z 12
    Create geotiff from existing map image around Barbados:
        python create_geotiff.py -s ~/Desktop/map_image.png --ULLR "13.1601374805, -59.6577740797, 13.1592965187, -59.6569159218"
""")
    parser.add_argument('--ULLR',
                        type=str,
                        help='Coordinates of upper-left and lower right coordinates of provided screenshot: "ULlat,ULlon,LRlat,LRlon"')
    parser.add_argument('--CLL',
                        type=str,
                        help='Center coordinates of google satellite imagery to be captured if we should create a new base image: "lat,lon"')
    parser.add_argument('-s', type=str, default='map_image.png', help='filepath of image to be turned into geotiff or location to store google satellite image')
    parser.add_argument('-z', type=str, default=20, help='zoom level of image to be downloaded from google, min zoom is 2 (whole earth), max zoom is 21')
    parser.add_argument('--pixel_size', default="640,640", type=str, help='Pixel size of image to be downloaded: "pixel_width,pixel_height" maximum size is 640x640')

    args = parser.parse_args()

    map_image = args.s
    zoom = args.z
    pixel_size_w, pixel_size_h = args.pixel_size.split(',')

    if args.CLL is not None:
        LC =  args.CLL.split(',')
        print("Received Center Coordinates, creating geotiff from google at %s" %map_image)
        UL, LR = download_screenshot(map_image, LC, pixel_size_w, pixel_size_h, zoom)
        create_geotiff(map_image, UL, LR)
    elif args.ULLR is not None:
        print("Received Corner Coordinates, creating geotiff from provided image: %s" %map_image)
        ULLR =  [x.strip() for x in args.ULLR.split(',')]
        UL = ULLR[:2]
        LR = ULLR[:-2]
        create_geotiff(map_image, UL, LR)