2

I have a point shapefile or feature class and I want to create bounding boxes around each point. The size of the bounding box will depend on user input, for example the user might want a bounding box of 50 x 50 meters.

The shapefiles are in WGS 84. This boils down to getting the x/y coordinates of a point x meters away from the existing point to the South East and to the North West (bearing 135 and 315 degrees).

How can I use arcpy to get the x/y coordinates of a point n meters away to the NW and the SE?

  • 2
    Why not just buffer the points by the desired amount and then use the Minimum Bounding Geometry tool to get the extent? – Evil Genius Aug 19 '15 at 13:58
  • A "box" in that coordinate reference shouldn't look rectangular (unless the point is on the equator), so you'll need at least four points (adding bearings 45 and 225). Large shapes (10+KM) should have even more vertices. This question is among the most common here at GIS SE – Vince Aug 19 '15 at 13:59
  • I am creating the box so I can query a wms service that requires a bounding box parameter. I don't expect the size to be more than 5 km in diameter, but most should be much smaller (50 to 500 m). – bdwyer Aug 19 '15 at 14:15
  • Based off Evil Genius' comment, I think I can run a buffer on the points and then get the arcpy geometry extent and use the XMin,XMax,YMin,YMax values to get my SE and NW points. – bdwyer Aug 19 '15 at 14:26
  • 1
    As long as the distances aren't massive and you don't need high levels of accuracy, you could just use trigonometry. You could even project the points on the fly if you wanted. – Paul Aug 19 '15 at 15:55
1

As @Paul suggested, you can use some simple trig to do this:

import arcpy

def getPoly(pt, w, h, wkid):
    """calculates polygon corners

    pt -- tuple of coordinates (x, y)
    w -- width of polygon (x length)
    h -- height of polygon (y length)
    sr -- spatial reference object
    """
    x, y = pt
    ll = (x - (w * 0.5), y - (h * 0.5))
    ul = (x - (w * 0.5), y + (h * 0.5))
    ur = (x + (w * 0.5), y + (h * 0.5))
    lr = (x + (w * 0.5), y - (h * 0.5))
    return arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in [ll,ul,ur,lr,ll]]), sr)

coords =(-10377711.530988107, 5617827.075599414)
sr = arcpy.SpatialReference(102100)
poly = getPoly(coords, 50, 50, sr) # make a bounding box 50x50 meters
e = poly.extent
bbox = ','.join(map(str, [e.XMin, e.YMin, e.XMax, e.YMax])) #here's the bbox

The bounding box came out as:

'-10377736.531,5617802.0756,-10377686.531,5617852.0756'

  • crmackey, thanks for the post. Will the code work with a shapefile in WGS84 (4326)? – bdwyer Aug 20 '15 at 15:12
  • Yep, it should work in any projection, just make sure you pass that WKID into the arcpy.SpatialReference() class. – crmackey Aug 20 '15 at 15:42

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.