diff --git a/salem/gis.py b/salem/gis.py
index d90a06a..a491af0 100644
--- a/salem/gis.py
+++ b/salem/gis.py
@@ -1223,13 +1223,18 @@ def transform_proj(p1, p2, x, y, nocopy=False):
         in case the two projections are equal, you can use nocopy if you wish
     """
 
-    if proj_is_same(p1, p2):
-        if nocopy:
-            return x, y
-        else:
-            return copy.deepcopy(x), copy.deepcopy(y)
+    try:
+        # This always makes a copy, even if projections are equivalent
+        return pyproj.transform(p1, p2, x, y,
+                                skip_equivalent=True, always_xy=True)
+    except TypeError:
+        if proj_is_same(p1, p2):
+            if nocopy:
+                return x, y
+            else:
+                return copy.deepcopy(x), copy.deepcopy(y)
 
-    return pyproj.transform(p1, p2, x, y)
+        return pyproj.transform(p1, p2, x, y)
 
 
 def transform_geometry(geom, crs=wgs84, to_crs=wgs84):
@@ -1470,7 +1475,7 @@ def mercator_grid(center_ll=None, extent=None, ny=600, nx=None,
     nx = np.rint(nx)
     ny = np.rint(ny)
 
-    e, n = pyproj.transform(wgs84, projloc, lon, lat)
+    e, n = transform_proj(wgs84, projloc, lon, lat)
 
     if origin== 'upper-left':
         corner = (-xx / 2. + e, yy / 2. + n)
@@ -1508,7 +1513,7 @@ def googlestatic_mercator_grid(center_ll=None, nx=640, ny=640, zoom=12, scale=1)
     xx = nx * mpix
     yy = ny * mpix
 
-    e, n = pyproj.transform(wgs84, projloc, lon, lat)
+    e, n = transform_proj(wgs84, projloc, lon, lat)
     corner = (-xx / 2. + e, yy / 2. + n)
     dxdy = (xx / nx, - yy / ny)
 
diff --git a/salem/tests/test_gis.py b/salem/tests/test_gis.py
index f874377..3f66741 100644
--- a/salem/tests/test_gis.py
+++ b/salem/tests/test_gis.py
@@ -967,50 +967,6 @@ def test_same_proj(self):
         if gis.has_gdal:
             self.assertTrue(gis.proj_is_same(p1, p2))
 
-    def test_pyproj_trafo(self):
-
-        x = np.random.randn(int(1e6)) * 60
-        y = np.random.randn(int(1e6)) * 60
-        t1 = time.time()
-        for i in np.arange(3):
-            xx, yy = pyproj.transform(wgs84, wgs84, x, y)
-        t1 = time.time() - t1
-        assert_allclose(xx, x)
-        assert_allclose(yy, y)
-
-        t2 = time.time()
-        for i in np.arange(3):
-            xx, yy = gis.transform_proj(wgs84, wgs84, x, y)
-        t2 = time.time() - t2
-        assert_allclose(xx, x)
-        assert_allclose(yy, y)
-
-        t3 = time.time()
-        for i in np.arange(3):
-            xx, yy = gis.transform_proj(wgs84, wgs84, x, y, nocopy=True)
-        t3 = time.time() - t3
-        assert_allclose(xx, x)
-        assert_allclose(yy, y)
-
-        self.assertTrue(t1 > t2)
-        self.assertTrue(t2 > t3)
-
-        t1 = time.time()
-        xx, yy = pyproj.transform(pyproj.Proj(init='epsg:26915'),
-                                  pyproj.Proj(init='epsg:26915'), x, y)
-        t1 = time.time() - t1
-        assert_allclose(xx, x, atol=1e-3)
-        assert_allclose(yy, y, atol=1e-3)
-
-        t2 = time.time()
-        xx, yy = gis.transform_proj(pyproj.Proj(init='epsg:26915'),
-                                    pyproj.Proj(init='epsg:26915'), x, y)
-        t2 = time.time() - t2
-        assert_allclose(xx, x)
-        assert_allclose(yy, y)
-
-        self.assertTrue(t1 > t2)
-
     @requires_shapely
     def test_geometry(self):
 
diff --git a/salem/tests/test_misc.py b/salem/tests/test_misc.py
index 69b25f9..6652bbe 100644
--- a/salem/tests/test_misc.py
+++ b/salem/tests/test_misc.py
@@ -197,6 +197,7 @@ def test_projplot(self):
 
         import pyproj
         import matplotlib.pyplot as plt
+        from salem.gis import transform_proj
 
         wgs84 = pyproj.Proj(proj='latlong', datum='WGS84')
         fig = plt.figure()
@@ -207,7 +208,7 @@ def test_projplot(self):
         proj_out = pyproj.Proj("+init=EPSG:4326", preserve_units=True)
         proj_in = pyproj.Proj(srs, preserve_units=True)
 
-        lon, lat = pyproj.transform(proj_in, proj_out, -2235000, -2235000)
+        lon, lat = transform_proj(proj_in, proj_out, -2235000, -2235000)
         np.testing.assert_allclose(lon, 70.75731, atol=1e-5)
 
 
diff --git a/salem/wrftools.py b/salem/wrftools.py
index 17bb6ee..81e2145 100644
--- a/salem/wrftools.py
+++ b/salem/wrftools.py
@@ -676,7 +676,7 @@ def geogrid_simulator(fpath, do_maps=True, map_kwargs=None):
     pwrf = gis.check_crs(pwrf)
 
     # get easting and northings from dom center (probably unnecessary here)
-    e, n = pyproj.transform(wgs84, pwrf, pargs['ref_lon'], pargs['lat_0'])
+    e, n = gis.transform_proj(wgs84, pwrf, pargs['ref_lon'], pargs['lat_0'])
 
     # LL corner
     nx, ny = e_we[0]-1, e_sn[0]-1