Getting a box

Sun 23 June 2013

Filed under Blog

Let's say one wants to find out a box around a given latitude/longitude point. How to do this without GIS software using an artisan UNIX approach? Hell, if it works for coffee! The steps are as follows:

Convert the centre point into UTM so we have meters rather than degrees
Write a box (manual buffering)
Convert the box corners back to WGS84

We need to know the projection we want to use as degrees are not constant at all latitudes. To know what UTM zone you need you could go here

For the example above of Lopé National Park, we could use UTM zone 32S with the WGS84 datum. This, according to spatialreference.org is EPSG code 32732. Let' now use gdaltransform to reproject to UTM:

$ echo 11.459789 -0.169737 | gdaltransform -s_srs "EPSG:4326" -t_srs "EPSG:32732"
773796.421559477 9981221.54175601 0

So this gives us the coordinates of the centre poin in meters. We know that we need a 10km square box around it, so we do this by defining the edges as (x-5000, y-5000), (x-5000, y+5000), (x+5000, y+5000), (x-5000, y+5000), (x-5000, y-5000). We can quickly just create a GML file with these points using a template like this:

<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
     xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
     xsi:schemaLocation="http://ogr.maptools.org/ gabon2.xsd"
     xmlns:ogr="http://ogr.maptools.org/"
     xmlns:gml="http://www.opengis.net/gml">

  <gml:featureMember>
    <ogr:Site fid="Site">
      <ogr:geometryProperty><gml:Polygon srsName="EPSG:32732"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>%g,%g,0 %g,%g,0 %g,%g,0 %g,%g,0 %g,%g,0</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
      <ogr:Name>Shape 1</ogr:Name>
      <ogr:Description></ogr:Description>
    </ogr:Site>
  </gml:featureMember>
</ogr:FeatureCollection>

In Python, we can just use the above snippet to populate the box, save it to a text file and convert it to our desired output:

gml = """<?xml version="1.0" encoding="utf-8" ?>
    <ogr:FeatureCollection
         xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
         xsi:schemaLocation="http://ogr.maptools.org/ gabon2.xsd"
         xmlns:ogr="http://ogr.maptools.org/"
         xmlns:gml="http://www.opengis.net/gml">

      <gml:featureMember>
        <ogr:Site fid="Site">
          <ogr:geometryProperty><gml:Polygon srsName="EPSG:32732"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>%g,%g,0 %g,%g,0 %g,%g,0 %g,%g,0 %g,%g,0</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
          <ogr:Name>Shape 1</ogr:Name>
          <ogr:Description></ogr:Description>
        </ogr:Site>
      </gml:featureMember>
    </ogr:FeatureCollection>"""
# The center location from gdaltransform
x,y = 773796.421559477, 9981221.54175601
fp = open ("my_box.gml", 'w' )
fp.write ( gml % ( x-5000, y-5000, x-5000, y+5000, x+5000, y+5000, x+5000, y-5000, x-5000, y-5000) )
fp.close()

Now, my_box.gml contains our box in UTM. What we want is a Shapefile in WGS84. We can use ogr2ogr for that:

$ ogr2ogr -f "ESRI Shapefile" -t_srs "EPSG:4326" my_box_wgs84.shp my_box.gml
$ ogr2ogr -f "KML" -t_srs "EPSG:4326" my_box_wgs84.kml my_box.gml

Check that we got something useful:

$ ogrinfo -al -so my_box_wgs84.kml
Had to open data source read-only.
INFO: Open of `my_box_wgs84.kml'
      using driver `KML' successful.

Layer name: Site
Geometry: Polygon
Feature Count: 1
Extent: (11.414887, -0.214953) - (11.504683, -0.124552)
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9108"]],
    AUTHORITY["EPSG","4326"]]
Name: String (0.0)
Description: String (0.0)

$ ogrinfo -al -so my_box_wgs84.shp
INFO: Open of `my_box_wgs84.shp'
      using driver `ESRI Shapefile' successful.

Layer name: my_box_wgs84
Geometry: Polygon
Feature Count: 1
Extent: (11.414887, -0.214953) - (11.504683, -0.124552)
Layer SRS WKT:
GEOGCS["GCS_WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS_84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]]
fid: String (80.0)
Name: String (7.0)

Getting a box

Comments