12

i have a shapefile of country subdivisions and i would like to extract an array of bounding lat and lon for each division.. is it possible to do?

3
  • Please clarify: do you have a shapefile or an Excel file?
    – whuber
    Commented Nov 24, 2010 at 15:08
  • 1
    The title doesn't really reflect your question, please think about editing it.
    – DavidF
    Commented Nov 24, 2010 at 15:21
  • i have a shapefile
    – mossplix
    Commented Nov 24, 2010 at 15:43

5 Answers 5

26

Using the ogr Python module from OSGEO, this example will give you a tuple containing the coords that define an envelope for each feature.

from osgeo import ogr

ds = ogr.Open("mn_counties.shp")
lyr = ds.GetLayerByName("mn_counties")

lyr.ResetReading()

for feat in lyr:
    # get bounding coords in minx, maxx, miny, maxy format
    env = feat.GetGeometryRef().GetEnvelope()
    # get bounding coords in minx, miny, maxx, maxy format
    bbox = [env[0], env[2], env[1], env[3]]
    print env
    print bbox
    print
1
  • 2
    ... and free geospatial python at that ; /
    – DavidF
    Commented Nov 24, 2010 at 15:56
4

One possible way to proceed using SAGA GIS http://www.saga-gis.org After opening your shapefile run these 3 modules: 1. Modules\Shapes\Tools\Get Shapes extent

  1. Modules\Shapes\Tools\Points\Points from lines [contrary to what the name suggest you can also use this to get points from a polygon]

  2. Modules\Shapes\Tools\Points\Add Coordinates to points This will give you a table which contains the x and y coordinates of the 4 corners of the bounding box of your polygon file.

4

In arcgis, here's python code. result is a list of minx, miny, maxx, maxy, minM, maxM, minZ, maxZ (

import arcpy
for feat in arcpy.SearchCursor(r"c:\data\f.gdb\counties"):
    print feat.Shape.extent

-2.66852727251546 49.4265363633626 -2.52848181818121 49.5079454546192 NaN NaN NaN NaN
-10.463336363782 51.4455454544593 -6.01305454583045 55.3799909091533 NaN NaN NaN NaN
-4.77778181827614 54.0555454544593 -4.35347272688468 54.4100000000002 NaN NaN NaN NaN
4

Here's an R version, using example data from the rgdal package:

library(rgdal)
dsn <- system.file("vectors/ps_cant_31.MIF", package = "rgdal")[1]
d <- readOGR(dsn = dsn, layer="ps_cant_31")

## transform if this is not longlat
if (is.projected(d)) d <- spTransform(d, CRS("+proj=longlat +ellps=WGS84"))

for (i in 1:nrow(d)) {
  print(bbox(d[i,]))    
}
1

I use fiona and shapely for that kind of tasks:

import fiona
from shapely.geometry import shape

with fiona.open(r'd:\Projects\_00_Data\_USstates\fe_2007_us_state00.shp', 'r') as features:
    for i, feat in enumerate(features):
        geom = shape(feat['geometry'])
        name = feat['properties']['NAME00']
        print ','.join((name,) + tuple([str(i) for i in geom.bounds]))
2
  • That does NOT provide Lat / Lon .
    – harpan
    Commented Dec 9, 2018 at 21:25
  • The output looks like this, those are the lat/lon coordinates: -124.72583900000001,45.544321,-116.915989,49.002494 -82.626182,37.202467,-77.71951899999999,40.638801 -111.056888,40.996345999999996,-104.052287,45.005904 -67.955811,17.913769,-65.22156799999999,18.511706
    – Matej
    Commented Jan 20, 2019 at 2:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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