How do I convert a GeoTiff's neatline metadata from vector format to pixels using GDAL?

Question

I have a GeoTIFF that produces the following when run through gdalinfo:

Driver: GTiff/GeoTIFF
Files: PAA_1.tif
Size is 14400, 6600
Coordinate System is:
PROJCS["PAA_1_Asia_South_Lambert_Conformal_Conic",
    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["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",-13.3333333],
    PARAMETER["standard_parallel_2",46.66666666],
    PARAMETER["latitude_of_origin",-15],
    PARAMETER["central_meridian",125],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]
GeoTransform =
  -7964169.743710292, 1177.758187333239, 322.1610457363416
  6367638.541192908, 322.2290109665665, -1177.712898399229
Metadata:
  AREA_OR_POINT=Area
  CREATION_DATE=D:20141215110023-05'00'
  CREATOR=ESRI ArcMap 9.3.1.4050
  NEATLINE=POLYGON ((-5337534.907769799232483 -155620.899077162524918,-5337535.365463367663324 10233860.720439316704869,10253410.032338818535209 10233861.175854712724686,10253410.596936786547303 -155621.599452315276721,-5337534.907769799232483 -155620.899077162524918))

My goal is to use the given neatline in the metadata to crop the image. I do not want to use gdalwarp because I do not want to rotate or warp the image, just crop it. Down the road, there may be other imaging tools (other than gdal) that crop the image, so I would prefer converting the given neatline coordinates into pixels relative to the top left corner and then pass those pixels off to something like imagemagick or PIL for image cropping.

Based on information from other answers, I put this together:

import gdal
from geometry import MapToPixel

ds = gdal.Open('PAA_1.tif')
transform = ds.GetGeoTransform()

# parses the neatline string into a list of (x,y) tuples
neatline = ds.GetMetadata()['NEATLINE'][10:-2]
neatline = [(pt.split(' ')[0], pt.split(' ')[1]) for pt in neatline.split(',')]
neatline = [(float(pt[0]), float(pt[1])) for pt in neatline]

# MapToPixel: https://code.google.com/p/metageta/source/browse/trunk/metageta/geometry.py#497
pixels = [MapToPixel(*pair, gt=transform) for pair in neatline]
print pixels

#RESULTS:
# [(665.3036244770069, 5720.95251596183), (2910.358041930786, -2486.5316409600555), (15226.421624842182, 883.213682625752), (12981.36704208365, 9090.698775704966), (665.3036244770069, 5720.95251596183)]

All this does is parse the NEATLINE metadata into x,y coordinate pairs and use the image's GeoTransform data, along with MapToPixel to convert the coordinates into pixel values. However, none of these values map to anything that looks like a neatline area of the image, and some of them are beyond the width and height bounds of the image.

If I understand correctly, it should be possible to apply some simple linear function to a set of coordinates to convert them into pixel values using the image's GeoTransform metadata. MapToPixel should be providing such a transformation. Please correct me if I'm wrong. I can also settle for hearing about GDAL utilities that could crop the image using NEATLINE without warping or rotating it (somehow preventing gdalwarp from doing this or using a different utility altogether).

Answer 1

I believe that you have done everything right but the neatline is crappy.

I checked the first coordinate of the neatline by using the formula for affine transformation from http://www.gdal.org/gdal_datamodel.html

Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2) 
Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)

The GeoTransform parameters of your image are

GT(0)= -7964169,744
GT(1)= 1177,758187
GT(2)= 322,1610457
GT(3)= 6367638,541
GT(4)= 322,229011
GT(5)= -1177,712898

First Xpixel, Yline pair from your formula is

(665.3036244770069, 5720.95251596183)

With these values Xgeo, Ygeo (rounded by Excel) is about the same as the first vertex in your neatline polygon:

-5337539.99, -155619.10 (Was -5337534.907769799232483, -155620.899077162524918)

I feat that you must make an intermediate copy of your map with gdalwarp and measure a better neatline. You are not the first one who has met neatline that does not cover the area of the map image in PDF gdalwarp not clipping neatline properly.

I do not know any GDAL utility that would support directly cropping by neatline but once you know the right pixels and the area to crop is rectangle without rotation you can use gdal_translate with -srswin

-srcwin xoff yoff xsize ysize: Selects a subwindow from the source image for copying based on pixel/line location.

Answer 2

You could create a polygon geometry directly from the neatline polygon WKT using ogr.CreateGeometryFromWkt(), then rasterize using gdal.RasterizeLayer().

# http://gis.stackexchange.com/a/131250/2856 (cc by-sa 3.0 attribution required)
import os
from osgeo import gdal,ogr
ogr.UseExceptions(True)

outfn='PPA_1_neatline.tif'

ds = gdal.Open('PAA_1.tif')
cols = ds.RasterXSize
rows = ds.RasterYSize
projection = ds.GetProjection()
geotransform = ds.GetGeoTransform()
bands = 1
datatype=gdal.GDT_Byte

neatline = ds.GetMetadata()['NEATLINE']

driver=ogr.GetDriverByName('Memory')
poly = ogr.CreateGeometryFromWkt(neatline)
memds=driver.CreateDataSource('tmpmemds')
lyr=memds.CreateLayer('neatline', geom_type=ogr.wkbPolygon)
feat = ogr.Feature(lyr.GetLayerDefn())
feat.SetGeometry(poly)
lyr.CreateFeature(feat)

driver = gdal.GetDriverByName('GTiff')
outds = driver.Create(outfn, cols, rows, bands, datatype)
outds.SetProjection(projection)
outds.SetGeoTransform(geotransform)
gdal.RasterizeLayer(outds, [1], lyr, None, None, [1], ['ALL_TOUCHED=TRUE'])