Merge pull request #11 from kjordahl/feature-transform

Add methods to transform coordinate systems

doc/docs.rst

@@ -208,6 +208,19 @@ Additionally, the following methods are implemented:
   Load a ``GeoSeries`` from a file from any format recognized by
   `fiona`_.
 
+.. method:: GeoSeries.to_crs(crs=None, epsg=None)
+
+  Transform all geometries in a GeoSeries to a different coordinate
+  reference system.  The ``crs`` attribute on the current GeoSeries
+  must be set.  Either ``crs`` in dictionary form or an EPSG code may
+  be specified for output.
+
+  This method will transform all points in all objects.  It has no
+  notion or projecting entire geometries.  All segments joining points
+  are assumed to be lines in the current projection, not geodesics.
+  Objects crossing the dateline (or other projection boundary) will
+  have undesirable behavior.
+
 .. method:: GeoSeries.plot(colormap='Set1')
 
   Generate a plot of the geometries in the ``GeoSeries``.
@@ -235,6 +248,21 @@ Currently only the following methods are implemented for a ``GeoDataFrame``:
   Load a ``GeoDataFrame`` from a file from any format recognized by
   `fiona`_.
 
+.. method:: GeoSeries.to_crs(crs=None, epsg=None, inplace=False)
+
+  Transform all geometries in the ``geometry`` column of a
+  GeoDataFrame to a different coordinate reference system.  The
+  ``crs`` attribute on the current GeoSeries must be set.  Either
+  ``crs`` in dictionary form or an EPSG code may be specified for
+  output.  If ``inplace=True`` the geometry column will be replaced in
+  the current dataframe, otherwise a new GeoDataFrame will be returned.
+
+  This method will transform all points in all objects.  It has no
+  notion or projecting entire geometries.  All segments joining points
+  are assumed to be lines in the current projection, not geodesics.
+  Objects crossing the dateline (or other projection boundary) will
+  have undesirable behavior.
+
 .. method:: GeoDataFrame.plot()
 
   Generate a plot of the geometries in the ``GeoDataFrame``.

geopandas/geodataframe.py

@@ -71,6 +71,25 @@ class GeoDataFrame(DataFrame):
              'features': [feature(i, row) for i, row in self.iterrows()]},
             **kwargs )
 
+    def to_crs(self, crs=None, epsg=None, inplace=False):
+        """Transform geometries to a new coordinate reference system
+
+        This method will transform all points in all objects.  It has
+        no notion or projecting entire geometries.  All segments
+        joining points are assumed to be lines in the current
+        projection, not geodesics.  Objects crossing the dateline (or
+        other projection boundary) will have undesirable behavior.
+        """
+        if inplace:
+            df = self
+        else:
+            df = self.copy()
+            df.crs = self.crs
+        geom = df.geometry.to_crs(crs=crs, epsg=epsg)
+        df.geometry = geom
+        if not inplace:
+            return df
+
     def __getitem__(self, key):
         """
         The geometry column is not stored as a GeoSeries, so need to make sure

geopandas/geoseries.py

@@ -1,13 +1,16 @@
 from warnings import warn
+from functools import partial
 
 import numpy as np
 from pandas import Series, DataFrame
 
+import pyproj
 from shapely.geometry import shape, Polygon, Point
 from shapely.geometry.collection import GeometryCollection
 from shapely.geometry.base import BaseGeometry
-from shapely.ops import cascaded_union, unary_union
+from shapely.ops import cascaded_union, unary_union, transform
 import fiona
+from fiona.crs import from_epsg
 
 from plotting import plot_series
 
@@ -380,3 +383,32 @@ class GeoSeries(Series):
 
     def plot(self, *args, **kwargs):
         return plot_series(self, *args, **kwargs)
+
+    #
+    # Additional methods
+    #
+
+    def to_crs(self, crs=None, epsg=None):
+        """Transform geometries to a new coordinate reference system
+
+        This method will transform all points in all objects.  It has
+        no notion or projecting entire geometries.  All segments
+        joining points are assumed to be lines in the current
+        projection, not geodesics.  Objects crossing the dateline (or
+        other projection boundary) will have undesirable behavior.
+        """
+        if self.crs is None:
+            raise ValueError('Cannot transform naive geometries.  '
+                             'Please set a crs on the object first.')
+        if crs is None:
+            try:
+                crs = from_epsg(epsg)
+            except TypeError:
+                raise TypeError('Must set either crs or epsg for output.')
+        proj_in = pyproj.Proj(**self.crs)
+        proj_out = pyproj.Proj(**crs)
+        project = partial(pyproj.transform, proj_in, proj_out)
+        result = self.apply(lambda geom: transform(project, geom))
+        result.__class__ = GeoSeries
+        result.crs = crs
+        return result

tests/test_geoseries.py

@@ -29,6 +29,10 @@ class TestSeries(unittest.TestCase):
         self.a1.index = ['A', 'B']
         self.a2 = self.g2.copy()
         self.a2.index = ['B', 'C']
+        self.esb = Point(-73.9847, 40.7484)
+        self.sol = Point(-74.0446, 40.6893)
+        self.landmarks = GeoSeries([self.esb, self.sol],
+                                   crs={'init': 'epsg:4326', 'no_defs': True})
 
     def test_area(self):
         assert_array_equal(self.g1.area.values, np.array([0.5, 1.0]))
@@ -177,3 +181,8 @@ class TestSeries(unittest.TestCase):
         self.assertTrue(np.alltrue(self.g2.contains(self.g2.representative_point())))
         self.assertTrue(np.alltrue(self.g3.contains(self.g3.representative_point())))
         self.assertTrue(np.alltrue(self.g4.contains(self.g4.representative_point())))
+
+    def test_transform(self):
+        utm18n = self.landmarks.to_crs(epsg=26918)
+        lonlat = utm18n.to_crs(epsg=4326)
+        self.assertTrue(np.alltrue(self.landmarks.almost_equals(lonlat)))