developmentseed/geolambda

createcopy causing coordinate system to change in geo tiff

Opened this issue · 3 comments

furick1 commented

I'm attempting to correct a compression issue with a geotiff by creating a copy or using translate which works locally with gdal_translate. The operation below using geolambda oddly causes the projected coordinate system to change from WGS 84 / Pseudo-Mercator to WGS_1984_Web_Mercator_Auxiliary_Sphere. I'm using the latest public geolambda layer referenced from arn:aws:lambda:us-east-1:552188055668:layer:geolambda:4 and geolambda-python layer from arn:aws:lambda:us-east-1:552188055668:layer:geolambda-python:3

Why is the PROJCS changing along with other values? Perhaps there is a better way to accomplish the task but still would like to know why the issue below happens.

Lambda Code

import logging
from osgeo import gdal
import boto3
from botocore.exceptions import ClientError

def upload_file(file_name, bucket, object_name=None):

# If S3 object_name was not specified, use file_name
if object_name is None:
    object_name = file_name

# Upload the file
s3_client = boto3.client('s3')
try:
    response = s3_client.upload_file(file_name, bucket, object_name)
except ClientError as e:
    logging.error(e)
    return False
return True

def lambda_handler(event, context=None):
    s3 = boto3.client('s3')

s3_filename = event["s3file"]
fname = event.get('filename', s3_filename)
fname = fname.replace('s3://', '/vsis3/')
bucket_name = 'bucket-name-here'
key = fname.replace('/vsis3/bucket-name-here/', '')
key = key.replace('.tif', '-gdal.tif')

dataset = gdal.Open(fname, gdal.GA_ReadOnly)
#gdal.Translate('/tmp/output.tif', dataset)

# Load the dataset into the virtual filesystem
temp_name = '/vsimem/current.tif'
tiff_driver = gdal.GetDriverByName('GTiff')
tiff_driver.CreateCopy(temp_name, dataset, False, [ 'COMPRESS=LZW' ])
#tiff_driver.Warp(temp_name,dataset,outputSRS='EPSG:4326')

# Read the raw data from the virtual filesystem
f = gdal.VSIFOpenL(temp_name, 'rb')
gdal.VSIFSeekL(f, 0, 2)  # seek to end
size = gdal.VSIFTellL(f)
gdal.VSIFSeekL(f, 0, 0)  # seek to beginning
data = gdal.VSIFReadL(1, size, f)
gdal.VSIFCloseL(f)

# Upload the raw data to s3
if s3.put_object(Key=key, Bucket=bucket_name, Body=data, ContentLength=size):
    details = 'GDAL conversion success. Input: ' + s3_filename + ' and Output: ' + key
    status = '200'
    summary_msg = 'success'
else:
    details = 'GDAL conversion failed. Input: ' + s3_filename + ' and Output: ' + key
    status = '400'
    summary_msg = 'error'

response = { 
    'status':status, 
    'message':{
        'summary' : summary_msg,
        'details' : details
    }
}
return response

gdal.Unlink(temp_name)
dataset = None
f = None

Input

{ "s3file": "s3://bucket-name/project/f228713b-cb88-41aa-8d0d-b728c3022d0e_2019:12:10T04:05:36:343_dsm.tif" }

Result

{ "status": "200", "message": { "summary": "success", "details": "GDAL conversion success. Input: s3://bucket-name/project/f228713b-cb88-41aa-8d0d-b728c3022d0e_2019:12:10T04:05:36:343_dsm.tif and Output: s3://bucket-name/project/f228713b-cb88-41aa-8d0d-b728c3022d0e_2019:12:10T04:05:36:343_dsm-gdal.tif" } }

Original gdalinfo

Driver: GTiff/GeoTIFF Files: f228713b-cb88-41aa-8d0d-b728c3022d0e_2019:12:10T04:05:36:343_dsm.tif Size is 15974, 14281 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, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], 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"]], AXIS["X",EAST], AXIS["Y",NORTH], 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 = (-8637199.014870001003146,4665140.998520000837743) Pixel Size = (0.020000000000000,-0.020000000000000) Metadata: AREA_OR_POINT=Area TIFFTAG_SOFTWARE=pix4dmapper Image Structure Metadata: COMPRESSION=LZW INTERLEAVE=BAND Corner Coordinates: Upper Left (-8637199.015, 4665140.999) ( 77d35'21.40"W, 38d36'15.29"N) Lower Left (-8637199.015, 4664855.379) ( 77d35'21.40"W, 38d36' 8.07"N) Upper Right (-8636879.535, 4665140.999) ( 77d35'11.07"W, 38d36'15.29"N) Lower Right (-8636879.535, 4664855.379) ( 77d35'11.07"W, 38d36' 8.07"N) Center (-8637039.275, 4664998.189) ( 77d35'16.24"W, 38d36'11.68"N) Band 1 Block=15974x1 Type=Float32, ColorInterp=Gray NoData Value=-10000

Output gdalinfo

Driver: GTiff/GeoTIFF Files: f228713b-cb88-41aa-8d0d-b728c3022d0e_2019:12:10T04:05:36:343_dsm-gdal.tif Size is 15974, 14281 Coordinate System is: PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere", 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 = (-8637199.014870001003146,4665140.998520000837743) Pixel Size = (0.020000000000000,-0.020000000000000) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left (-8637199.015, 4665140.999) ( 77d35'21.40"W, 38d36'15.29"N) Lower Left (-8637199.015, 4664855.379) ( 77d35'21.40"W, 38d36' 8.07"N) Upper Right (-8636879.535, 4665140.999) ( 77d35'11.07"W, 38d36'15.29"N) Lower Right (-8636879.535, 4664855.379) ( 77d35'11.07"W, 38d36' 8.07"N) Center (-8637039.275, 4664998.189) ( 77d35'16.24"W, 38d36'11.68"N) Band 1 Block=15974x1 Type=Float32, ColorInterp=Gray

gdal_translate local command that works

gdal_translate -q f228713b-cb88-41aa-8d0d-b728c3022d0e_2019:12:10T04:05:36:343_dsm.tif f228713b-cb88-41aa-8d0d-b728c3022d0e_2019:12:10T04:05:36:343_dsm-gdal.tif -co COMPRESS=LZW

furick1 commented

I just tried the older versions of the public geolambda layers, 2 and 1 respectively, and the same code works and doesn't change the coordinate system. hmmm

matthewhanson commented

@furick1 I'm not sure, but I wonder if it's due to the use of PROJ.6 now in GDAL 3 rather than PROJ.4 which was a huge refactor in how it worked.

I'm not entirely clear on the differences between Web Mercator and Web Mercator Aux Sphere, but it looks like there really isn't much difference at all. I wonder if PROJ always used WMAS but didn't necessarily label it as such, and so this was a clarification.

Have you found at any more about this?

furick1 commented