2

I'm using GDAL in Python to reproject rasters. I'd like to be able to take an existing raster, and project it into a new raster to my projection and output pixel size of choice.

I'm worried about how to calculate the bounding box in projected coordinates that I will in turn base the geotransform off of. gdal has a function called gdal.AutoCreateWarpedVRT that apparently creates a memory raster that I can then directly read the geotransform from. But I sometimes have VERY large rasters that don't even fit in main memory, I'd hate to create a memory based raster just to get the geotransform.

Is it sufficient to transform the vertices of the original bounding box into the projected coordinate system and then somehow make that box a square? (I'm sure the box will be distorted.)

|improve this question
  • 2
    Yes (with conditions). Take the bounding box, regardless of what coordinate system the raster currently is in it will be a rectangle, then project the 4 coordinates to give your output extent... sort of. Some projections swell in the centre (like geographic to mercator) in these cases extra points are needed in the middle of the box; to overcome odd shapes at least add a point on the edge in the middle of each of the corners though it would be best to densify (add a vertex every 10 pixels or so) the lines then get the envelope min/max from all the points. – Michael Stimson Oct 1 '15 at 22:16
  • Thanks @MichaelMiles-Stimson, I ended up doing what you suggested; making a 8 point bounding box (4 on the corners, for in the middle of the edges), reprojecting each point, then defining the bounding box as a box that captures all 8 transformed points. Any chance you can post this answer so I can mark it correct? – Rich Oct 9 '15 at 20:27
  • Pretty much what AndreJ said in his answer. It would be nice to have some working code though, Andre is using command line utilities which can be done using subprocess.popen; regardless of method the result is the same. At the moment I'm too busy write, code and reference an answer so I would support this one or write your own with code example. – Michael Stimson Oct 11 '15 at 22:23
  • Thanks @michael-miles-stimson, I posted the code I ended up writing that works on any number of edge sample points. – Rich Jun 6 '16 at 21:59
  • I guess this depends on what you're using the bounding box for, but if you're checking whether points lie within it, then it might make more sense to make a function that projects the points into the CRS, then checks them against the orthogonal box. – naught101 Mar 13 at 6:11
3

A safe way to get the exact new extent is to

  1. create a polygon from the raster extent with gdaltindex
  2. densify the polygon geometry
  3. reproject the polygon to the target CRS with ogr2ogr
  4. calculate the extent in the new CRS
|improve this answer
  • Thanks @AndreJ, if you don't mind, I posted what I ended up using which procedurally samples each edge of the bounding polygon, projects those points, then finds the tightest bounding box there. It's a little more data light since it doesn't need to densify a polygon nor make any intermediate files and also works directly in Python/GDAL. – Rich Jun 6 '16 at 22:07
1

It's been a few months since I posted this. @michael-miles-stimson suggested that the correct way to do this would be to sample points along each edge of the bounding box, project them, then build an encompassing bounding box on those points. The drop-in Python/GDAL code below does this and allows you to select the number of interpolating points you'd like to sample along the edge.

import numpy
from osgeo import osr


def transform_bounding_box(
        bounding_box, base_epsg, new_epsg, edge_samples=11):
    """Transform input bounding box to output projection.

    This transform accounts for the fact that the reprojected square bounding
    box might be warped in the new coordinate system.  To account for this,
    the function samples points along the original bounding box edges and
    attempts to make the largest bounding box around any transformed point
    on the edge whether corners or warped edges.

    Parameters:
        bounding_box (list): a list of 4 coordinates in `base_epsg` coordinate
            system describing the bound in the order [xmin, ymin, xmax, ymax]
        base_epsg (int): the EPSG code of the input coordinate system
        new_epsg (int): the EPSG code of the desired output coordinate system
        edge_samples (int): the number of interpolated points along each
            bounding box edge to sample along. A value of 2 will sample just
            the corners while a value of 3 will also sample the corners and
            the midpoint.

    Returns:
        A list of the form [xmin, ymin, xmax, ymax] that describes the largest
        fitting bounding box around the original warped bounding box in
        `new_epsg` coordinate system.
    """
    base_ref = osr.SpatialReference()
    base_ref.ImportFromEPSG(base_epsg)

    new_ref = osr.SpatialReference()
    new_ref.ImportFromEPSG(new_epsg)

    transformer = osr.CoordinateTransformation(base_ref, new_ref)

    p_0 = numpy.array((bounding_box[0], bounding_box[3]))
    p_1 = numpy.array((bounding_box[0], bounding_box[1]))
    p_2 = numpy.array((bounding_box[2], bounding_box[1]))
    p_3 = numpy.array((bounding_box[2], bounding_box[3]))

    def _transform_point(point):
        trans_x, trans_y, _ = (transformer.TransformPoint(*point))
        return (trans_x, trans_y)

    # This list comprehension iterates over each edge of the bounding box,
    # divides each edge into `edge_samples` number of points, then reduces
    # that list to an appropriate `bounding_fn` given the edge.
    # For example the left edge needs to be the minimum x coordinate so
    # we generate `edge_samples` number of points between the upper left and
    # lower left point, transform them all to the new coordinate system
    # then get the minimum x coordinate "min(p[0] ...)" of the batch.
    transformed_bounding_box = [
        bounding_fn(
            [_transform_point(
                p_a * v + p_b * (1 - v)) for v in numpy.linspace(
                    0, 1, edge_samples)])
        for p_a, p_b, bounding_fn in [
            (p_0, p_1, lambda p_list: min([p[0] for p in p_list])),
            (p_1, p_2, lambda p_list: min([p[1] for p in p_list])),
            (p_2, p_3, lambda p_list: max([p[0] for p in p_list])),
            (p_3, p_0, lambda p_list: max([p[1] for p in p_list]))]]
    return transformed_bounding_box
|improve this answer

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.