2

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).

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.

|improve this answer|||||
  • Many thanks for the cross checking. I wondered if the source of data was bad, but it's good to have some confidence. I will load up the image in QGIS as mentioned in the answer you referenced for neatline checking. Your gdal_translate method of cropping with a good image and the pixel values in-hand is helpful too. Thank you! – Sam Atcheson Jan 21 '15 at 16:13
0

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'])
|improve this answer|||||
  • Should memds=driver.CreateDataSource(outshape) be memds=driver.CreateDataSource(outfn) instead? outshape is undefined. When running this code with the outfn fix, I was getting black TIFs on output, but it's likely a bad cutline issue as mentioned in the other answer. I'll need to try this out on some other image to see if it works, but it seems like a good idea. – Sam Atcheson Jan 22 '15 at 21:30
  • @SamAtcheson No, outfn is the output tif. The polygon dataset is stored in memory. I had an outshape variable as I originally tested using a shapefile, then changed to the in memory driver. I must've deleted the variable definition when changing working code so it matched your data. Answer edited to correct. – user2856 Jan 23 '15 at 1:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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