Repository: geopandas/geopandas

Issue:
ENH: move CRS metadata to array level (allow columns with different CRS)
Currently, the CRS (`.crs` attribute) is kept track of at the GeoDataFrame (or GeoSeries if selecting the column) level. This means it is implicitly the CRS for the "active" geometry column. 

We could also look into moving the metadata to the GeometryArray. This feels more correct anyhow, and would allow to have multiple columns with different CRS (which would be especially relevant if returning a GeoSeries for other geometry columns as well -> https://github.com/geopandas/geopandas/issues/1192). 
The GeoDataFrame.crs would keep pointing to the crs of the active geometry column, to ensure to preserve the current behaviour.

Discussion:
martinfleis: I've fixed tests to make them compatible with older pandas, but there is one failing with pandas master, which I am not sure if should be fixed here or in pandas (I am not sure which PR causes it).

```py
s1 = GeoSeries(
    [
        Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]),
        Polygon([(2, 2), (4, 2), (4, 4), (2, 4)]),
    ]
)
df1 = GeoDataFrame({"col1": [1, 2], "geometry": s1}, crs="EPSG:4326")
```

This should **not** alter CRS of `s1` and `s1.values`, but with pandas master it does set CRS to `s1.values` during creation of `df1`. It looks that until few days ago, pandas made a copy before assigning array to df, but now it does not. @jorisvandenbossche any ideas?
martinfleis: I believe that DEV CI issue is pandas related. Opened an issue there: https://github.com/pandas-dev/pandas/issues/32960
martinfleis: Thanks! It is a bit messy at the moment and there might be bits and pieces which are no longer needed. 

> all operations that return new geometries (like centroid, buffer, intersection, etc) should preserve the CRS

Good point, I'll make sure this works on the array level. 
martinfleis: This seems to be working now.

Summary:
Moved CRS to GeometryArray level. Each column can then have its own CRS, not only one for whole GeoDataFrame. By default, `GeoSeries.crs` and `GeoDataFrame.crs` are linked directly to underlying GeometryArray (active one). In case of gdf, it also stores it as metadata, because we still allow GeoDataFrame without geometry column. In case of following workflow, crs has to be stored somewhere in the meantime:

```py
df = GeoDataFrame(crs=4326)
df = df.set_geometry(array)
```

One thing can be a bit tricky in terms of backward compatibility. That is the priority of CRS. Currently, if GeoSeries passed as geometry have CRS, gdf's CRS is ignored. I have mirrored this logic so current priority is `GeometryArray > GeoSeries > GeoDataFrame`. The closer we are to the geometry the higher probability is that the CRS is correct and should be kept. Examples:

```py
geoms = [Point(0, 0), Point(1, 1)]
array = from_shapely(geoms, crs=4326)
array.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

series = GeoSeries(array, crs=27700) # 27700 is ignored
series.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

df = GeoDataFrame(geometry=series, crs=5355)  # 5355 is ignored
df.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
```

Are you okay with this behaviour? Should we raise a warning in case of difference between set and passed CRS?
jorisvandenbossche: Thanks a lot for this effort!

> In case of gdf, it also stores it as metadata, because we still allow GeoDataFrame without geometry column. In case of following workflow, crs has to be stored somewhere in the meantime:
>
> ```
> df = GeoDataFrame(crs=4326)
>df = df.set_geometry(array)
> ```

Generally I don't like very much the "geodataframe without geometry column" case, and this is not making it better ;) I am wondering if we can't deprecate this: not speaking about the general aspect here (that's for another issue), but just the about passing the crs in GeoDataFrame if there are no geometry values yet. 
As for the above use case, the user can also do `df = GeoDataFrame(..); df.set_geometry(array, crs=4326)`, right?

We still need to store the crs at the GeoDataFrame level for now, but only for this specific use case.

> One thing can be a bit tricky in terms of backward compatibility. That is the priority of CRS. 

The rule of "The closer we are to the geometry the higher probability is that the CRS is correct and should be kept." sounds sensible. But on the other hand, if we forget for a moment the current behaviour, and look at `GeoSeries(array, crs=27700)` or `GeoDataFrame(geometry=series, crs=5355)`, what would you expect? I think I would expect that passing this `crs=` in the GeoSeries or GeoDataFrame constructor would *override* the crs of the geometries (not knowing the current behaviour). 
But, that's maybe also a good reason to at least warn about it

---

There might also be some other corner cases, like what to do when "overriding" existing geometries via assignment. Eg assume you already have a valid GeoDataFrame with a "geometry" column, and then do:

```
gdf['geometry'] = array-like
```

which is not uncommon when doing eg a map or apply, basically anything where you don't use a vectorized operation but loop over the geometry scalars. 

Should that try to "preserve" the crs from `gdf['geometry']` in case if the array-like has no crs information? 
From pandas point of view, when doing `gdf['col']` you create a whole new column, that has no connection with the previous one (eg it will also not try to preserve the dtype if that would be possible, like `df['existing_float_col'] = 1` will an int column). 
So the most logical thing might be to follow pandas in this, but that might also be a backwards compatible break if the crs is no longer stored "globally".


martinfleis: @jorisvandenbossche 

> I don't like very much the "geodataframe without geometry column" case

I don't like it either and I would vote for deprecation of such an option altogether (return pandas DataFrame in that case). I will add deprecation warning here that we will not support assigning CRS to GeoDataFrame without geometry column. We may also open an issue (unless there is one, I remember discussion on this) to discuss how to deprecate GeoDataFrame without geometry altogether.

> As for the above use case, the user can also do `df = GeoDataFrame(..); df.set_geometry(array, crs=4326)`, right?

Yes, that is preferable option as CRS is directly linked with geometry here.

> look at `GeoSeries(array, crs=27700)` or `GeoDataFrame(geometry=series, crs=5355)`, what would you expect?

I guess that without thinking about it, I would expect that it gets overwritten. But I don't think it is a behaviour we should go for. GeoPandas should be aware of CRS even if user is not, so in this case I would ideally warn that there is a mismatch and one should use `to_crs()`. 

We are also a bit inconsistent here, because `GeoDataFrame(geometry=series, crs=5355)` does not overwrite it, but `df.set_geometry(series, crs=5355)` does (it is clearly stated in docstring though). I would expect that this will be essentially the same thing, so the same behaviour. 

> Should that try to "preserve" the crs from `gdf['geometry']` in case if the array-like has no crs information?

This is a tricky case. It would be logical to follow pandas case, but keep in mind that our `'geometry'` column has a specific position. It is not just a column. If I want to, for example. keep only first 3 coordinates of LineString using `gdf['geometry'] = gdf.geometry.apply(...)` I would certainly not expect that I will lose CRS during the operation. So yes, I would try to preserve CRS even though it may not be fully consistent with pandas, because it is a common workflow here.
jorisvandenbossche: > This is a tricky case. It would be logical to follow pandas case, but keep in mind that our 'geometry' column has a specific position. It is not just a column. If I want to, for example, keep only first 3 coordinates of LineString using `gdf['geometry'] = gdf.geometry.apply(...)` I would certainly not expect that I will lose CRS during the operation. So yes, I would try to preserve CRS even though it may not be fully consistent with pandas, because it is a common workflow here.

We can try to have some checks for this, but we need to be aware that there are many corner cases to this:

- What if you assign it to a new column `gdf['other_geometry] = df.geometry.apply(...)` ? We could in principle still check for that.
- What if another geometry column was used like `gdf['other_geometry].apply(..)`, while this other column can have another CRS now than the overall dataframe one?
- What if you actually transformed the the geometries in the apply?

Now, many of those problems we also have today, of course (although the second actually becomes a problem), and is not specific to this PR. But it's a good moment to consider those.

We could also put some logic in `apply` to check if the output are still geometries and if so preserve the crs (potentially with an option to not do this). That could remediate part of it for when using `apply`. But of course not when doing a manual loop like `df['geometry'] = [.. for geom in df.geometry]`


martinfleis: I like the idea of _clever_ `apply`, adapting it to preserve CRS with the ability to opt-out/assign different CRS. Manual loop is a bit more tricky, as resulting array currently inherits GeoDataFrame's CRS (there is no other option). We can give it df's CRS, but that would not work if it is a result of a loop over another column. None of the solutions seems to be optimal.

I would go with following:

1. adapt `apply` to handle CRS as another PR
2. assign gdf's CRS to loop-based array only if it is passed to a `'geometry'` column (that should be the current behaviour of this PR) as I am not sure if we actually have another option
3. add a check of CRS mismatch between different levels of hierarchy (`GeometryArray > GeoSeries > GeoDataFrame`) within this PR
4. implement `.set_crs()` as more straightforward method of overwriting CRS (as another PR) mirroring `to_crs()` structure.
5. properly document all these cases (initially within this PR, updated by following PRs).

Does it make sense?
jorisvandenbossche: That does all makes sense to me, yes.

> assign gdf's CRS to loop-based array only if it is passed to a 'geometry' column (that should be the current behaviour of this PR) as I am not sure if we actually have another option

For loop-based assignment to columns that didn't yet exist / were not yet geometry columns, we right now also don't make it geometry-dtype columns (it will be object columns), so for those the question is even irrelevant (unless we would start inferring in `GeoDataFrame.__setitem__` if the values are potentially all shapely scalars, but I don't think we want to get into this)
martinfleis: I have added few commits preserving CRS for loop over active `'geometry'` column to preserve current behaviour and to check CRS match during creation of GeoSeries and GeoDataFrame. Plus short section of docs.

I have also opened two new issues as a follow-up. Ideally should be implemented before 0.8 release to have it all together.

@jorisvandenbossche could you review the code once you'll find a moment so we can get closer to merge? Thank you!
jorisvandenbossche: Regarding the "priority of the CRS", I was experimenting a bit more with master, and it seems we are actually not that consistent. Eg with GeoDataframe, the crs keyword gets ignored in some cases, but not all, and GeoSeries also does not do this.

```python
In [1]: import geopandas 
In [2]: import shapely.geometry 
In [3]: geoseries = geopandas.GeoSeries([shapely.geometry.Point(1, 1)], crs="EPSG:4326")

In [5]: geoseries.crs
Out[5]: 
<Geographic 2D CRS: EPSG:4326>
...

# passing a GeoSeries to GeoSeries(..) constructor -> the crs actually gets ignored ...
# (which is a bug I would say, but this might already be fixed on this PR?)
In [6]: geopandas.GeoSeries(geoseries).crs

# so specifying a CRS in the GeoSeries constructor overrules the passed geometries
In [8]: geopandas.GeoSeries(geoseries, crs="EPSG:31370").crs  
Out[8]: 
<Projected CRS: EPSG:31370>
...

# GeoDataFrame constructor -> passing geometries to "geometry" keyword
# -> the crs keyword gets ignored, and the result uses CRS of geometries
In [9]: geopandas.GeoDataFrame({'data': [1]}, geometry=geoseries, crs="EPSG:31370").crs
Out[9]: 
<Geographic 2D CRS: EPSG:4326>
...

# however, when passing the geomeries inside the dictionary, now the CRS of the
# geometries is ignored, and the "crs" keyword is used ..
In [10]: geopandas.GeoDataFrame({'data': [1], 'geometry': geoseries}, crs="EPSG:31370").crs 
Out[10]: 
<Projected CRS: EPSG:31370>
...

# and when passing inside the dictionary without specifying the crs keyword
# the crs gets lost (similar to series)
In [13]: geopandas.GeoDataFrame({'data': [1], 'geometry': geoseries}).crs 

# creating a geodataframe with EPSG:4326
In [14]: gdf = geopandas.GeoDataFrame({'data': [1]}, geometry=geoseries)  

In [15]: gdf.crs 
Out[15]: 
<Geographic 2D CRS: EPSG:4326>
...

# passing a gdf through GeoDataFrame constructor -> crs gets lost
In [18]: geopandas.GeoDataFrame(gdf).crs

# specifying the crs keyword -> that gets used
In [19]: geopandas.GeoDataFrame(gdf, crs="EPSG:31370").crs  
Out[19]: 
<Projected CRS: EPSG:31370>
...
```

So it all seems rather inconsistent ... There is indeed the case where `GeoDataFrame(geometry=.., crs=..)` where `crs` keyword gets ignored. But so that's certainly not universally the case.


jorisvandenbossche: Thanks for the updates and answers! One quick reply: 

> > I am still not really convinced this is the best option
>
> What would you prefer instead? Overwrite CRS? Automatically re-project? I believe that this is the safest option. You don't want to accidentally overwrite existing CRS

Yeah, it's certainly true that overwriting an existing CRS without warning is also not good .. 
But so if actually don't want users to do this, we should maybe make the warning into an error? (or start with a warning but say this will raise an error in the future) Then we force the user to be explicit in setting it afterwards?


martinfleis: > we should maybe make the warning into an error?

That is a good idea. That way we make sure that user is aware what is going on with the CRS and that it did not get overwritten. I'll change into FutureWarning now and we can switch to error in the (not so distant) future.
martinfleis: I just noticed that test of CRS during `__setitem__` started failing with pandas master recently. @jorisvandenbossche do you have any idea what causes it and on which side it should be fixed? Maybe https://github.com/pandas-dev/pandas/pull/33273?. On April 2 it was okay, April 6 is already failing. (https://travis-ci.org/github/geopandas/geopandas/jobs/671808342#L3570)
jorisvandenbossche: Hmm, not that PR I think, but will take a look.
jorisvandenbossche: So it might be the same PR as was mentioned in https://github.com/pandas-dev/pandas/issues/32960, but that is already from longer ago ..
martinfleis: You may be right. I have added a patch for pandas master to make sure that CRS is not lost during `__setitem__`. 
jorisvandenbossche: Related, and this might conflict or fix (not sure :-)) with the problem on pandas master, but we need to check the CRS in GeometryArray's `__setitem__`. If the value that is being set (eg for a slice) is a GeometryArray, we should check that the CRS matches (if one is set). 

See also https://github.com/pandas-dev/pandas/issues/33457, which I suppose is the reason for the failure.
martinfleis: >  we need to check the CRS in GeometryArray's `__setitem__`

If we check in GA and assign it there, it fixes the issue with pandas master. I just have to make sure that all options are covered and test it. Nice catch.
jorisvandenbossche: > If we check in GA and assign it there, it fixes the issue with pandas master. 

Yeah, I only think that in principle, we should *not* assign the crs then, since you can only set a part of the geometries (or at least we should only allow to override the crs in ``_setitem__`` with eg a ``slice(None)`` from `array[:] = ...` so you know you are replacing all geometries). 
But since that is a bit tricky, I think in principle we should error. But that will only be possible once the regression is fixed on pandas' side.
martinfleis: >  I think in principle we should error. But that will only be possible once the regression is fixed on pandas' side.

I agree with raising error. I have prepared check (commented out) for that and reverted previous fix as it may have caused incorrect behaviour. Now CRS during `__setitem__` is failing with pandas master and I would wait for pandas to fix that. There are also other consequences of that regression, like following, which raises CRS mismatch even though both are the same:

```py
gs = GeoSeries(polys, crs=27700)
GeoDataFrame({"geometry": gs, "df1": [0, 1, 2]}, crs=27700)

FutureWarning: CRS mismatch between CRS of the passed geometries and 'crs'. Use 'GeoDataFrame.crs = crs' to overwrite CRS or 'GeoDataFrame.to_crs()' to reproject geometries. CRS mismatch will raise an error in the future versions of GeoPandas.
  GeoDataFrame({"geometry": gs, "df1": [0, 1, 2]}, crs=27700)
```
Something just get lost during some operations in pandas master.
jorisvandenbossche: Thanks a lot for the answers and updates! 

I think all remaining items have follow-up issues, so let's get this in?
martinfleis: Couple of minor updates per discussion in.

> I think all remaining items have follow-up issues, so let's get this in?

Yes, I have opened issues for all mentioned items. We can move on.
martinfleis: One more thing. The `__setitem__` issue with pandas master discussed above. Now it is failing as we loose CRS during pandas operations. Are we okay with it as it is likely pandas issue? I may open an issue for that as well to keep an eye on it.
jorisvandenbossche: Yes, let's create an issue for that as well

Patch:
diff --git a/doc/source/projections.rst b/doc/source/projections.rst
index d6febeec..9903c041 100644
--- a/doc/source/projections.rst
+++ b/doc/source/projections.rst
@@ -98,6 +98,44 @@ Re-projecting is the process of changing the representation of locations from on
     ax.set_title("Mercator");
 
 
+Projection for multiple geometry columns
+----------------------------------------
+
+GeoPandas 0.8 implements support for different projections assigned to different geometry
+columns of the same GeoDataFrame. The projection is now stored together with geometries per column (directly
+on the GeometryArray level).
+
+Note that if GeometryArray has assigned projection, it is preferred over the
+projection passed to GeoSeries or GeoDataFrame during the creation:
+
+.. code-block:: python
+
+   >>> array.crs
+   <Geographic 2D CRS: EPSG:4326>
+   Name: WGS 84
+   Axis Info [ellipsoidal]:
+   - Lat[north]: Geodetic latitude (degree)
+   - Lon[east]: Geodetic longitude (degree)
+   ...
+   >>> GeoSeries(array, crs=3395).crs  # crs=3395 is ignored as array already has CRS
+   FutureWarning: CRS mismatch between CRS of the passed geometries and 'crs'. Use 'GeoSeries.crs = crs' to overwrite CRS or 'GeoSeries.to_crs()' to reproject geometries. CRS mismatch will raise an error in the future versions of GeoPandas.
+       GeoSeries(array, crs=3395).crs
+
+   <Geographic 2D CRS: EPSG:4326>
+   Name: WGS 84
+   Axis Info [ellipsoidal]:
+   - Lat[north]: Geodetic latitude (degree)
+   - Lon[east]: Geodetic longitude (degree)
+   ...
+
+If you want to overwrite projection, you can then assign it to the GeoSeries manually
+or re-project geometries to the target projection using either ``GeoSeries.crs = 3395``
+or ``GeoSeries.to_crs(3395)``.
+
+All GeometryArray-based operations preserve projection; however, if you loop over a column
+containing geometry, this information might be lost.
+
+
 Upgrading to GeoPandas 0.7 with pyproj > 2.2 and PROJ > 6
 ---------------------------------------------------------
 
@@ -384,7 +422,7 @@ If we construct the CRS object from the EPSG code (truncated output):
 
 You can see that the CRS object constructed from the WKT string has a "Easting,
 Northing" (i.e. x, y) axis order, while the CRS object constructed from the EPSG
-code has a (Northing, Easting) axis order. 
+code has a (Northing, Easting) axis order.
 
 Only having this difference in axis order is no problem when using the CRS in
 GeoPandas, since GeoPandas always uses a (x, y) order to store the data
diff --git a/geopandas/array.py b/geopandas/array.py
index 21deaf03..44c3546b 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -12,6 +12,7 @@ import shapely.geometry
 from shapely.geometry.base import BaseGeometry
 import shapely.ops
 import shapely.wkt
+from pyproj import CRS
 
 try:
     import pygeos
@@ -98,13 +99,23 @@ def _is_scalar_geometry(geom):
         return isinstance(geom, BaseGeometry)
 
 
-def from_shapely(data):
+def from_shapely(data, crs=None):
     """
     Convert a list or array of shapely objects to a GeometryArray.
 
     Validates the elements.
+
+    Parameters
+    ----------
+    data : array-like
+        list or array of shapely objects
+    crs : value, optional
+        Coordinate Reference System of the geometry objects. Can be anything accepted by
+        :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+        such as an authority string (eg "EPSG:4326") or a WKT string.
+
     """
-    return GeometryArray(vectorized.from_shapely(data))
+    return GeometryArray(vectorized.from_shapely(data), crs=crs)
 
 
 def to_shapely(geoms):
@@ -116,11 +127,21 @@ def to_shapely(geoms):
     return vectorized.to_shapely(geoms.data)
 
 
-def from_wkb(data):
+def from_wkb(data, crs=None):
     """
     Convert a list or array of WKB objects to a GeometryArray.
+
+    Parameters
+    ----------
+    data : array-like
+        list or array of WKB objects
+    crs : value, optional
+        Coordinate Reference System of the geometry objects. Can be anything accepted by
+        :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+        such as an authority string (eg "EPSG:4326") or a WKT string.
+
     """
-    return GeometryArray(vectorized.from_wkb(data))
+    return GeometryArray(vectorized.from_wkb(data), crs=crs)
 
 
 def to_wkb(geoms):
@@ -132,11 +153,21 @@ def to_wkb(geoms):
     return vectorized.to_wkb(geoms.data)
 
 
-def from_wkt(data):
+def from_wkt(data, crs=None):
     """
     Convert a list or array of WKT objects to a GeometryArray.
+
+    Parameters
+    ----------
+    data : array-like
+        list or array of WKT objects
+    crs : value, optional
+        Coordinate Reference System of the geometry objects. Can be anything accepted by
+        :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+        such as an authority string (eg "EPSG:4326") or a WKT string.
+
     """
-    return GeometryArray(vectorized.from_wkt(data))
+    return GeometryArray(vectorized.from_wkt(data), crs=crs)
 
 
 def to_wkt(geoms, **kwargs):
@@ -148,13 +179,17 @@ def to_wkt(geoms, **kwargs):
     return vectorized.to_wkt(geoms.data, **kwargs)
 
 
-def points_from_xy(x, y, z=None):
+def points_from_xy(x, y, z=None, crs=None):
     """
     Generate GeometryArray of shapely Point geometries from x, y(, z) coordinates.
 
     Parameters
     ----------
     x, y, z : iterable
+    crs : value, optional
+        Coordinate Reference System of the geometry objects. Can be anything accepted by
+        :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+        such as an authority string (eg "EPSG:4326") or a WKT string.
 
     Examples
     --------
@@ -167,7 +202,7 @@ def points_from_xy(x, y, z=None):
     -------
     output : GeometryArray
     """
-    return GeometryArray(vectorized.points_from_xy(x, y, z))
+    return GeometryArray(vectorized.points_from_xy(x, y, z), crs=crs)
 
 
 class GeometryArray(ExtensionArray):
@@ -178,7 +213,7 @@ class GeometryArray(ExtensionArray):
 
     _dtype = GeometryDtype()
 
-    def __init__(self, data):
+    def __init__(self, data, crs=None):
         if isinstance(data, self.__class__):
             data = data.data
         elif not isinstance(data, np.ndarray):
@@ -192,6 +227,28 @@ class GeometryArray(ExtensionArray):
             )
         self.data = data
 
+        self._crs = None
+        self.crs = crs
+
+    @property
+    def crs(self):
+        """
+        The Coordinate Reference System (CRS) represented as a ``pyproj.CRS``
+        object.
+
+        Returns None if the CRS is not set, and to set the value it
+        :getter: Returns a ``pyproj.CRS`` or None. When setting, the value
+        Coordinate Reference System of the geometry objects. Can be anything accepted by
+        :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+        such as an authority string (eg "EPSG:4326") or a WKT string.
+        """
+        return self._crs
+
+    @crs.setter
+    def crs(self, value):
+        """Sets the value of the crs"""
+        self._crs = None if not value else CRS.from_user_input(value)
+
     @property
     def dtype(self):
         return self._dtype
@@ -208,7 +265,7 @@ class GeometryArray(ExtensionArray):
             # to numpy array, pass-through non-array-like indexers
             idx = pd.api.indexers.check_array_indexer(self, idx)
         if isinstance(idx, (Iterable, slice)):
-            return GeometryArray(self.data[idx])
+            return GeometryArray(self.data[idx], crs=self.crs)
         else:
             raise TypeError("Index type not supported", idx)
 
@@ -245,13 +302,22 @@ class GeometryArray(ExtensionArray):
                 "Value should be either a BaseGeometry or None, got %s" % str(value)
             )
 
+        # TODO: use this once pandas-dev/pandas#33457 is fixed
+        # if hasattr(value, "crs"):
+        #     if value.crs and (value.crs != self.crs):
+        #         raise ValueError(
+        #             "CRS mismatch between CRS of the passed geometries "
+        #             "and CRS of existing geometries."
+        #         )
+
     if compat.USE_PYGEOS:
 
         def __getstate__(self):
-            return pygeos.to_wkb(self.data)
+            return (pygeos.to_wkb(self.data), self._crs)
 
         def __setstate__(self, state):
-            geoms = pygeos.from_wkb(state)
+            geoms = pygeos.from_wkb(state[0])
+            self._crs = state[1]
             self.data = geoms
             self.base = None
 
@@ -301,23 +367,23 @@ class GeometryArray(ExtensionArray):
 
     @property
     def boundary(self):
-        return GeometryArray(vectorized.boundary(self.data))
+        return GeometryArray(vectorized.boundary(self.data), crs=self.crs)
 
     @property
     def centroid(self):
-        return GeometryArray(vectorized.centroid(self.data))
+        return GeometryArray(vectorized.centroid(self.data), crs=self.crs)
 
     @property
     def convex_hull(self):
-        return GeometryArray(vectorized.convex_hull(self.data))
+        return GeometryArray(vectorized.convex_hull(self.data), crs=self.crs)
 
     @property
     def envelope(self):
-        return GeometryArray(vectorized.envelope(self.data))
+        return GeometryArray(vectorized.envelope(self.data), crs=self.crs)
 
     @property
     def exterior(self):
-        return GeometryArray(vectorized.exterior(self.data))
+        return GeometryArray(vectorized.exterior(self.data), crs=self.crs)
 
     @property
     def interiors(self):
@@ -325,7 +391,7 @@ class GeometryArray(ExtensionArray):
         return vectorized.interiors(self.data)
 
     def representative_point(self):
-        return GeometryArray(vectorized.representative_point(self.data))
+        return GeometryArray(vectorized.representative_point(self.data), crs=self.crs)
 
     #
     # Binary predicates
@@ -382,16 +448,22 @@ class GeometryArray(ExtensionArray):
     #
 
     def difference(self, other):
-        return GeometryArray(self._binary_method("difference", self, other))
+        return GeometryArray(
+            self._binary_method("difference", self, other), crs=self.crs
+        )
 
     def intersection(self, other):
-        return GeometryArray(self._binary_method("intersection", self, other))
+        return GeometryArray(
+            self._binary_method("intersection", self, other), crs=self.crs
+        )
 
     def symmetric_difference(self, other):
-        return GeometryArray(self._binary_method("symmetric_difference", self, other))
+        return GeometryArray(
+            self._binary_method("symmetric_difference", self, other), crs=self.crs
+        )
 
     def union(self, other):
-        return GeometryArray(self._binary_method("union", self, other))
+        return GeometryArray(self._binary_method("union", self, other), crs=self.crs)
 
     #
     # Other operations
@@ -402,19 +474,22 @@ class GeometryArray(ExtensionArray):
 
     def buffer(self, distance, resolution=16, **kwargs):
         return GeometryArray(
-            vectorized.buffer(self.data, distance, resolution=resolution, **kwargs)
+            vectorized.buffer(self.data, distance, resolution=resolution, **kwargs),
+            crs=self.crs,
         )
 
     def interpolate(self, distance, normalized=False):
         return GeometryArray(
-            vectorized.interpolate(self.data, distance, normalized=normalized)
+            vectorized.interpolate(self.data, distance, normalized=normalized),
+            crs=self.crs,
         )
 
     def simplify(self, tolerance, preserve_topology=True):
         return GeometryArray(
             vectorized.simplify(
                 self.data, tolerance, preserve_topology=preserve_topology
-            )
+            ),
+            crs=self.crs,
         )
 
     def project(self, other, normalized=False):
@@ -442,33 +517,38 @@ class GeometryArray(ExtensionArray):
 
     def affine_transform(self, matrix):
         return GeometryArray(
-            vectorized._affinity_method("affine_transform", self.data, matrix)
+            vectorized._affinity_method("affine_transform", self.data, matrix),
+            crs=self.crs,
         )
 
     def translate(self, xoff=0.0, yoff=0.0, zoff=0.0):
         return GeometryArray(
-            vectorized._affinity_method("translate", self.data, xoff, yoff, zoff)
+            vectorized._affinity_method("translate", self.data, xoff, yoff, zoff),
+            crs=self.crs,
         )
 
     def rotate(self, angle, origin="center", use_radians=False):
         return GeometryArray(
             vectorized._affinity_method(
                 "rotate", self.data, angle, origin=origin, use_radians=use_radians
-            )
+            ),
+            crs=self.crs,
         )
 
     def scale(self, xfact=1.0, yfact=1.0, zfact=1.0, origin="center"):
         return GeometryArray(
             vectorized._affinity_method(
                 "scale", self.data, xfact, yfact, zfact, origin=origin
-            )
+            ),
+            crs=self.crs,
         )
 
     def skew(self, xs=0.0, ys=0.0, origin="center", use_radians=False):
         return GeometryArray(
             vectorized._affinity_method(
                 "skew", self.data, xs, ys, origin=origin, use_radians=use_radians
-            )
+            ),
+            crs=self.crs,
         )
 
     #
@@ -531,7 +611,7 @@ class GeometryArray(ExtensionArray):
 
     def copy(self, *args, **kwargs):
         # still taking args/kwargs for compat with pandas 0.24
-        return GeometryArray(self.data.copy())
+        return GeometryArray(self.data.copy(), crs=self._crs)
 
     def take(self, indices, allow_fill=False, fill_value=None):
         from pandas.api.extensions import take
@@ -547,7 +627,7 @@ class GeometryArray(ExtensionArray):
         result = take(self.data, indices, allow_fill=allow_fill, fill_value=fill_value)
         if allow_fill and fill_value is None:
             result[pd.isna(result)] = None
-        return GeometryArray(result)
+        return GeometryArray(result, crs=self.crs)
 
     def _fill(self, idx, value):
         """ Fill index locations with value
@@ -802,7 +882,7 @@ class GeometryArray(ExtensionArray):
         ExtensionArray
         """
         data = np.concatenate([ga.data for ga in to_concat])
-        return GeometryArray(data)
+        return GeometryArray(data, crs=to_concat[0].crs)
 
     def _reduce(self, name, skipna=True, **kwargs):
         # including the base class version here (that raises by default)
diff --git a/geopandas/base.py b/geopandas/base.py
index 610e310f..9305f3cb 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -4,7 +4,6 @@ import numpy as np
 import pandas as pd
 from pandas import DataFrame, MultiIndex, Series
 
-from pyproj import CRS
 from shapely.geometry import box
 from shapely.geometry.base import BaseGeometry
 from shapely.ops import cascaded_union
@@ -146,12 +145,12 @@ class GeoPandasBase(object):
         :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
         such as an authority string (eg "EPSG:4326") or a WKT string.
         """
-        return self._crs
+        return self.geometry.values.crs
 
     @crs.setter
     def crs(self, value):
         """Sets the value of the crs"""
-        self._crs = None if not value else CRS.from_user_input(value)
+        self.geometry.values.crs = value
 
     @property
     def geom_type(self):
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 617d98c6..c5c3ec8d 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -1,4 +1,5 @@
 import json
+import warnings
 
 import numpy as np
 import pandas as pd
@@ -6,6 +7,8 @@ from pandas import DataFrame, Series
 
 from shapely.geometry import mapping, shape
 
+from pyproj import CRS
+
 from geopandas.array import GeometryArray, from_shapely, GeometryDtype
 from geopandas.base import GeoPandasBase, is_geometry_type
 from geopandas.geoseries import GeoSeries
@@ -15,12 +18,14 @@ from geopandas.plotting import plot_dataframe
 DEFAULT_GEO_COLUMN_NAME = "geometry"
 
 
-def _ensure_geometry(data):
+def _ensure_geometry(data, crs=None):
     """
     Ensure the data is of geometry dtype or converted to it.
 
     If input is a (Geo)Series, output is a GeoSeries, otherwise output
     is GeometryArray.
+
+    If the input is a GeometryDtype with a set CRS, `crs` is ignored.
     """
     if is_geometry_type(data):
         if isinstance(data, Series):
@@ -28,10 +33,10 @@ def _ensure_geometry(data):
         return data
     else:
         if isinstance(data, Series):
-            out = from_shapely(np.asarray(data))
+            out = from_shapely(np.asarray(data), crs=crs)
             return GeoSeries(out, index=data.index, name=data.name)
         else:
-            out = from_shapely(data)
+            out = from_shapely(data, crs=crs)
             return out
 
 
@@ -43,8 +48,10 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
 
     Parameters
     ----------
-    crs : str (optional)
-        Coordinate system
+    crs : value (optional)
+        Coordinate Reference System of the geometry objects. Can be anything accepted by
+        :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+        such as an authority string (eg "EPSG:4326") or a WKT string.
     geometry : str or array (optional)
         If str, column to use as geometry. If array, will be set as 'geometry'
         column on GeoDataFrame.
@@ -61,8 +68,7 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
 
         # need to set this before calling self['geometry'], because
         # getitem accesses crs
-        self._crs = None
-        self.crs = crs
+        self._crs = crs if crs is not None else None
 
         # set_geometry ensures the geometry data have the proper dtype,
         # but is not called if `geometry=None` ('geometry' column present
@@ -74,7 +80,23 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
             # only if we have actual geometry values -> call set_geometry
             index = self.index
             try:
-                self["geometry"] = _ensure_geometry(self["geometry"].values)
+                if (
+                    hasattr(self["geometry"].values, "crs")
+                    and self["geometry"].values.crs
+                    and crs
+                    and not self["geometry"].values.crs == crs
+                ):
+                    warnings.warn(
+                        "CRS mismatch between CRS of the passed geometries "
+                        "and 'crs'. Use 'GeoDataFrame.crs = crs' to overwrite CRS "
+                        "or 'GeoDataFrame.to_crs()' to reproject geometries. "
+                        "CRS mismatch will raise an error in the future versions "
+                        "of GeoPandas.",
+                        FutureWarning,
+                        stacklevel=2,
+                    )  # TODO: change 'GeoDataFrame.crs = crs' to 'set_crs()' once done
+                    # TODO: raise error in 0.9 or 0.10.
+                self["geometry"] = _ensure_geometry(self["geometry"].values, crs)
             except TypeError:
                 pass
             else:
@@ -86,9 +108,33 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
                 geometry = "geometry"
 
         if geometry is not None:
+            if (
+                hasattr(geometry, "crs")
+                and geometry.crs
+                and crs
+                and not geometry.crs == crs
+            ):
+                warnings.warn(
+                    "CRS mismatch between CRS of the passed geometries "
+                    "and 'crs'. Use 'GeoDataFrame.crs = crs' to overwrite CRS "
+                    "or 'GeoDataFrame.to_crs()' to reproject geometries. "
+                    "CRS mismatch will raise an error in the future versions "
+                    "of GeoPandas.",
+                    FutureWarning,
+                    stacklevel=2,
+                )  # TODO: change 'GeoDataFrame.crs = crs' to 'set_crs()' once done
+                # TODO: raise error in 0.9 or 0.10.
             self.set_geometry(geometry, inplace=True)
         self._invalidate_sindex()
 
+        if geometry is None and crs:
+            warnings.warn(
+                "Assigning CRS to a GeoDataFrame without a geometry column is now "
+                "deprecated and will not be supported in the future.",
+                FutureWarning,
+                stacklevel=2,
+            )
+
     def __setattr__(self, attr, val):
         # have to special case geometry b/c pandas tries to use as column...
         if attr == "geometry":
@@ -149,9 +195,6 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
         else:
             frame = self.copy()
 
-        if not crs:
-            crs = getattr(col, "crs", self.crs)
-
         to_remove = None
         geo_column_name = self._geometry_column_name
         if isinstance(col, (Series, list, np.ndarray, GeometryArray)):
@@ -174,13 +217,17 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
         if to_remove:
             del frame[to_remove]
 
-        if isinstance(level, GeoSeries) and level.crs != crs:
+        if not crs:
+            level_crs = getattr(level, "crs", None)
+            crs = level_crs if level_crs is not None else self._crs
+
+        if isinstance(level, (GeoSeries, GeometryArray)) and level.crs != crs:
             # Avoids caching issues/crs sharing issues
             level = level.copy()
             level.crs = crs
 
         # Check that we are using a listlike of geometries
-        level = _ensure_geometry(level)
+        level = _ensure_geometry(level, crs=crs)
         index = frame.index
         frame[geo_column_name] = level
         if frame.index is not index and len(frame.index) == len(index):
@@ -225,6 +272,39 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
         self.rename(columns={geometry_col: col}, inplace=inplace)
         self.set_geometry(col, inplace=inplace)
 
+    @property
+    def crs(self):
+        """
+        The Coordinate Reference System (CRS) represented as a ``pyproj.CRS``
+        object.
+
+        Returns None if the CRS is not set, and to set the value it
+        :getter: Returns a ``pyproj.CRS`` or None. When setting, the value
+        can be anything accepted by
+        :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+        such as an authority string (eg "EPSG:4326") or a WKT string.
+        """
+        return self._crs
+
+    @crs.setter
+    def crs(self, value):
+        """Sets the value of the crs"""
+        if self._geometry_column_name not in self:
+            warnings.warn(
+                "Assigning CRS to a GeoDataFrame without a geometry column is now "
+                "deprecated and will not be supported in the future.",
+                FutureWarning,
+                stacklevel=4,
+            )
+            self._crs = None if not value else CRS.from_user_input(value)
+        else:
+            if hasattr(self.geometry.values, "crs"):
+                self.geometry.values.crs = value
+                self._crs = self.geometry.values.crs
+            else:
+                # column called 'geometry' without geometry
+                self._crs = None if not value else CRS.from_user_input(value)
+
     @classmethod
     def from_file(cls, filename, **kwargs):
         """Alternate constructor to create a ``GeoDataFrame`` from a file.
@@ -543,13 +623,13 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
     def to_crs(self, crs=None, epsg=None, inplace=False):
         """Transform geometries to a new coordinate reference system.
 
-        Transform all geometries in a GeoSeries to a different coordinate
+        Transform all geometries in an active geometry column to a different coordinate
         reference system.  The ``crs`` attribute on the current GeoSeries must
         be set.  Either ``crs`` or ``epsg`` may be specified for output.
 
-        This method will transform all points in all objects.  It has no notion
+        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
+        assumed to be lines in the current projection, not geodesics. Objects
         crossing the dateline (or other projection boundary) will have
         undesirable behavior.
 
@@ -589,17 +669,27 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
         geo_col = self._geometry_column_name
         if isinstance(result, Series) and isinstance(result.dtype, GeometryDtype):
             result.__class__ = GeoSeries
-            result.crs = self.crs
             result._invalidate_sindex()
         elif isinstance(result, DataFrame) and geo_col in result:
             result.__class__ = GeoDataFrame
-            result.crs = self.crs
             result._geometry_column_name = geo_col
             result._invalidate_sindex()
         elif isinstance(result, DataFrame) and geo_col not in result:
             result.__class__ = DataFrame
         return result
 
+    def __setitem__(self, key, value):
+        """
+        Overwritten to preserve CRS of GeometryArray in cases like
+        df['geometry'] = [geom... for geom in df.geometry]
+        """
+        if not pd.api.types.is_list_like(key) and key == self._geometry_column_name:
+            try:
+                value = _ensure_geometry(value, crs=self.crs)
+            except TypeError:
+                warnings.warn("Geometry column does not contain geometry.")
+        super(GeoDataFrame, self).__setitem__(key, value)
+
     #
     # Implement pandas methods
     #
@@ -653,6 +743,7 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
         else:
             for name in self._metadata:
                 object.__setattr__(self, name, getattr(other, name, None))
+
         return self
 
     def plot(self, *args, **kwargs):
diff --git a/geopandas/geoseries.py b/geopandas/geoseries.py
index 44afb1ba..e0f194a1 100644
--- a/geopandas/geoseries.py
+++ b/geopandas/geoseries.py
@@ -52,8 +52,11 @@ class GeoSeries(GeoPandasBase, Series):
         The geometries to store in the GeoSeries.
     index : array-like or Index
         The index for the GeoSeries.
-    crs : str, dict (optional)
-        Coordinate Reference System of the geometry objects.
+    crs : value (optional)
+        Coordinate Reference System of the geometry objects. Can be anything accepted by
+        :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+        such as an authority string (eg "EPSG:4326") or a WKT string.
+
     kwargs
         Additional arguments passed to the Series constructor,
          e.g. ``name``.
@@ -76,11 +79,29 @@ class GeoSeries(GeoPandasBase, Series):
 
     """
 
-    _metadata = ["name", "_crs"]
+    _metadata = ["name"]
 
     def __new__(cls, data=None, index=None, crs=None, **kwargs):
         # we need to use __new__ because we want to return Series instance
         # instead of GeoSeries instance in case of non-geometry data
+
+        if hasattr(data, "crs") and crs:
+            if not data.crs:
+                # make a copy to avoid setting CRS to passed GeometryArray
+                data = data.copy()
+            else:
+                if not data.crs == crs:
+                    warnings.warn(
+                        "CRS mismatch between CRS of the passed geometries "
+                        "and 'crs'. Use 'GeoSeries.crs = crs' to overwrite CRS "
+                        "or 'GeoSeries.to_crs()' to reproject geometries. "
+                        "CRS mismatch will raise an error in the future versions "
+                        "of GeoPandas.",
+                        FutureWarning,
+                        stacklevel=2,
+                    )  # TODO: change 'GeoSeries.crs = crs' to 'set_crs()' once done
+                    # TODO: raise error in 0.9 or 0.10.
+
         if isinstance(data, SingleBlockManager):
             if isinstance(data.blocks[0].dtype, GeometryDtype):
                 if data.blocks[0].ndim == 2:
@@ -96,7 +117,7 @@ class GeoSeries(GeoPandasBase, Series):
                     data = SingleBlockManager([block], data.axes[0], fastpath=True)
                 self = super(GeoSeries, cls).__new__(cls)
                 super(GeoSeries, self).__init__(data, index=index, **kwargs)
-                self.crs = crs
+                self.crs = getattr(self.values, "crs", crs)
                 return self
             warnings.warn(_SERIES_WARNING_MSG, FutureWarning, stacklevel=2)
             return Series(data, index=index, **kwargs)
@@ -125,7 +146,7 @@ class GeoSeries(GeoPandasBase, Series):
                     return s
             # try to convert to GeometryArray, if fails return plain Series
             try:
-                data = from_shapely(s.values)
+                data = from_shapely(s.values, crs)
             except TypeError:
                 warnings.warn(_SERIES_WARNING_MSG, FutureWarning, stacklevel=2)
                 return s
@@ -134,7 +155,9 @@ class GeoSeries(GeoPandasBase, Series):
 
         self = super(GeoSeries, cls).__new__(cls)
         super(GeoSeries, self).__init__(data, index=index, name=name, **kwargs)
-        self.crs = crs
+
+        if not self.crs:
+            self.crs = crs
         self._invalidate_sindex()
         return self
 
diff --git a/geopandas/tools/overlay.py b/geopandas/tools/overlay.py
index 33ec74be..8d0b7c43 100644
--- a/geopandas/tools/overlay.py
+++ b/geopandas/tools/overlay.py
@@ -86,7 +86,7 @@ def _overlay_difference(df1, df2):
             lambda x, y: x.difference(y), [geom] + list(df2.geometry.iloc[neighbours])
         )
         new_g.append(new)
-    differences = GeoSeries(new_g, index=df1.index)
+    differences = GeoSeries(new_g, index=df1.index, crs=df1.crs)
     poly_ix = differences.type.isin(["Polygon", "MultiPolygon"])
     differences.loc[poly_ix] = differences[poly_ix].buffer(0)
     geom_diff = differences[~differences.is_empty].copy()


Test patch:
diff --git a/geopandas/testing.py b/geopandas/testing.py
index 5c20b13e..901059b7 100644
--- a/geopandas/testing.py
+++ b/geopandas/testing.py
@@ -212,14 +212,16 @@ def assert_geodataframe_equal(
     )
 
     # geometry comparison
-    assert_geoseries_equal(
-        left.geometry,
-        right.geometry,
-        check_dtype=check_dtype,
-        check_less_precise=check_less_precise,
-        check_geom_type=check_geom_type,
-        check_crs=False,
-    )
+    for col, dtype in left.dtypes.iteritems():
+        if isinstance(dtype, GeometryDtype):
+            assert_geoseries_equal(
+                left[col],
+                right[col],
+                check_dtype=check_dtype,
+                check_less_precise=check_less_precise,
+                check_geom_type=check_geom_type,
+                check_crs=check_crs,
+            )
 
     # drop geometries and check remaining columns
     left2 = left.drop([left._geometry_column_name], axis=1)
diff --git a/geopandas/tests/test_crs.py b/geopandas/tests/test_crs.py
index a0667d93..088cc37b 100644
--- a/geopandas/tests/test_crs.py
+++ b/geopandas/tests/test_crs.py
@@ -1,12 +1,16 @@
 from distutils.version import LooseVersion
 import os
 
+import random
+
 import numpy as np
+import pandas as pd
 
-from shapely.geometry import Point
+from shapely.geometry import Point, Polygon
 import pyproj
 
-from geopandas import GeoDataFrame, points_from_xy
+from geopandas import GeoSeries, GeoDataFrame, points_from_xy, datasets, read_file
+from geopandas.array import from_shapely, from_wkb, from_wkt
 
 from geopandas.testing import assert_geodataframe_equal
 import pytest
@@ -133,3 +137,399 @@ def test_skip_exact_same():
     df = df_epsg26918()
     utm = df.to_crs(df.crs)
     assert_geodataframe_equal(df, utm, check_less_precise=True)
+
+
+# Test CRS on GeometryArray level
+class TestGeometryArrayCRS:
+    def setup_method(self):
+        self.osgb = pyproj.CRS(27700)
+        self.wgs = pyproj.CRS(4326)
+
+        self.geoms = [Point(0, 0), Point(1, 1)]
+        self.polys = [
+            Polygon([(random.random(), random.random()) for i in range(3)])
+            for _ in range(10)
+        ]
+        self.arr = from_shapely(self.polys, crs=27700)
+
+    def test_array(self):
+        arr = from_shapely(self.geoms)
+        arr.crs = 27700
+        assert arr.crs == self.osgb
+
+        arr = from_shapely(self.geoms, crs=27700)
+        assert arr.crs == self.osgb
+
+    def test_series(self):
+        s = GeoSeries(crs=27700)
+        assert s.crs == self.osgb
+        assert s.values.crs == self.osgb
+
+        arr = from_shapely(self.geoms)
+        s = GeoSeries(arr, crs=27700)
+        assert s.crs == self.osgb
+        assert s.values.crs == self.osgb
+
+        # manually change CRS
+        s.crs = 4326
+        assert s.crs == self.wgs
+        assert s.values.crs == self.wgs
+
+        s = GeoSeries(self.geoms, crs=27700)
+        assert s.crs == self.osgb
+        assert s.values.crs == self.osgb
+
+        arr = from_shapely(self.geoms, crs=27700)
+        s = GeoSeries(arr)
+        assert s.crs == self.osgb
+        assert s.values.crs == self.osgb
+
+        with pytest.warns(FutureWarning):
+            s = GeoSeries(arr, crs=4326)
+        assert s.crs == self.osgb
+
+    @pytest.mark.filterwarnings("ignore:Assigning CRS")
+    def test_dataframe(self):
+        arr = from_shapely(self.geoms, crs=27700)
+        df = GeoDataFrame(geometry=arr)
+        assert df.crs == self.osgb
+        assert df.geometry.crs == self.osgb
+        assert df.geometry.values.crs == self.osgb
+
+        arr = from_shapely(self.geoms)
+        s = GeoSeries(arr, crs=27700)
+        df = GeoDataFrame(geometry=s)
+        assert df.crs == self.osgb
+        assert df.geometry.crs == self.osgb
+        assert df.geometry.values.crs == self.osgb
+
+        # different passed CRS than array CRS is ignored
+        with pytest.warns(FutureWarning):
+            df = GeoDataFrame(geometry=s, crs=4326)
+        assert df.crs == self.osgb
+        assert df.geometry.crs == self.osgb
+        assert df.geometry.values.crs == self.osgb
+        with pytest.warns(FutureWarning):
+            GeoDataFrame(geometry=s, crs=4326)
+        with pytest.warns(FutureWarning):
+            GeoDataFrame({"data": [1, 2], "geometry": s}, crs=4326)
+        with pytest.warns(FutureWarning):
+            GeoDataFrame(df, crs=4326).crs
+
+        # manually change CRS
+        arr = from_shapely(self.geoms)
+        s = GeoSeries(arr, crs=27700)
+        df = GeoDataFrame(geometry=s)
+        df.crs = 4326
+        assert df.crs == self.wgs
+        assert df.geometry.crs == self.wgs
+        assert df.geometry.values.crs == self.wgs
+
+        df = GeoDataFrame(self.geoms, columns=["geom"], crs=27700)
+        assert df.crs == self.osgb
+        df = df.set_geometry("geom")
+        assert df.crs == self.osgb
+        assert df.geometry.crs == self.osgb
+        assert df.geometry.values.crs == self.osgb
+        assert df.geom.crs == self.osgb
+        assert df.geom.values.crs == self.osgb
+
+        df = GeoDataFrame(geometry=self.geoms, crs=27700)
+        assert df.crs == self.osgb
+        assert df.geometry.crs == self.osgb
+        assert df.geometry.values.crs == self.osgb
+
+        df = GeoDataFrame(crs=27700)
+        df = df.set_geometry(self.geoms)
+        assert df.crs == self.osgb
+        assert df.geometry.crs == self.osgb
+        assert df.geometry.values.crs == self.osgb
+
+        # new geometry with set CRS has priority over GDF CRS
+        df = GeoDataFrame(crs=27700)
+        df = df.set_geometry(self.geoms, crs=4326)
+        assert df.crs == self.wgs
+        assert df.geometry.crs == self.wgs
+        assert df.geometry.values.crs == self.wgs
+
+        df = GeoDataFrame()
+        df = df.set_geometry(s)
+        assert df.crs == self.osgb
+        assert df.geometry.crs == self.osgb
+        assert df.geometry.values.crs == self.osgb
+
+        arr = from_shapely(self.geoms, crs=27700)
+        df = GeoDataFrame()
+        df = df.set_geometry(arr)
+        assert df.crs == self.osgb
+        assert df.geometry.crs == self.osgb
+        assert df.geometry.values.crs == self.osgb
+
+        arr = from_shapely(self.geoms)
+        df = GeoDataFrame({"col1": [1, 2], "geometry": arr}, crs=4326)
+        assert df.crs == self.wgs
+        assert df.geometry.crs == self.wgs
+        assert df.geometry.values.crs == self.wgs
+
+        arr = from_shapely(self.geoms, crs=4326)
+        df = GeoDataFrame({"col1": [1, 2], "geometry": arr})
+        assert df.crs == self.wgs
+        assert df.geometry.crs == self.wgs
+        assert df.geometry.values.crs == self.wgs
+
+        # geometry column without geometry
+        df = GeoDataFrame({"geometry": [0, 1]})
+        df.crs = 27700
+        assert df.crs == self.osgb
+
+    def test_read_file(self):
+        nybb_filename = datasets.get_path("nybb")
+        df = read_file(nybb_filename)
+        assert df.crs == pyproj.CRS(2263)
+        assert df.geometry.crs == pyproj.CRS(2263)
+        assert df.geometry.values.crs == pyproj.CRS(2263)
+
+    def test_multiple_geoms(self):
+        arr = from_shapely(self.geoms, crs=27700)
+        s = GeoSeries(self.geoms, crs=4326)
+        df = GeoDataFrame(s, geometry=arr, columns=["col1"])
+        assert df.crs == self.osgb
+        assert df.geometry.crs == self.osgb
+        assert df.geometry.values.crs == self.osgb
+        assert df.col1.crs == self.wgs
+        assert df.col1.values.crs == self.wgs
+
+    def test_multiple_geoms_set_geom(self):
+        arr = from_shapely(self.geoms, crs=27700)
+        s = GeoSeries(self.geoms, crs=4326)
+        df = GeoDataFrame(s, geometry=arr, columns=["col1"])
+        df = df.set_geometry("col1")
+        assert df.crs == self.wgs
+        assert df.geometry.crs == self.wgs
+        assert df.geometry.values.crs == self.wgs
+        assert df["geometry"].crs == self.osgb
+        assert df["geometry"].values.crs == self.osgb
+
+    def test_assign_cols(self):
+        arr = from_shapely(self.geoms, crs=27700)
+        s = GeoSeries(self.geoms, crs=4326)
+        df = GeoDataFrame(s, geometry=arr, columns=["col1"])
+        df["geom2"] = s
+        df["geom3"] = s.values
+        df["geom4"] = from_shapely(self.geoms)
+        assert df.crs == self.osgb
+        assert df.geometry.crs == self.osgb
+        assert df.geometry.values.crs == self.osgb
+        assert df.geom2.crs == self.wgs
+        assert df.geom2.values.crs == self.wgs
+        assert df.geom3.crs == self.wgs
+        assert df.geom3.values.crs == self.wgs
+        assert df.geom4.crs is None
+        assert df.geom4.values.crs is None
+
+    def test_copy(self):
+        arr = from_shapely(self.geoms, crs=27700)
+        s = GeoSeries(self.geoms, crs=4326)
+        df = GeoDataFrame(s, geometry=arr, columns=["col1"])
+
+        arr_copy = arr.copy()
+        assert arr_copy.crs == arr.crs
+
+        s_copy = s.copy()
+        assert s_copy.crs == s.crs
+        assert s_copy.values.crs == s.values.crs
+
+        df_copy = df.copy()
+        assert df_copy.crs == df.crs
+        assert df_copy.geometry.crs == df.geometry.crs
+        assert df_copy.geometry.values.crs == df.geometry.values.crs
+        assert df_copy.col1.crs == df.col1.crs
+        assert df_copy.col1.values.crs == df.col1.values.crs
+
+    def test_rename(self):
+        arr = from_shapely(self.geoms, crs=27700)
+        s = GeoSeries(self.geoms, crs=4326)
+        df = GeoDataFrame(s, geometry=arr, columns=["col1"])
+        df = df.rename(columns={"geometry": "geom"}).set_geometry("geom")
+        assert df.crs == self.osgb
+        assert df.geometry.crs == self.osgb
+        assert df.geometry.values.crs == self.osgb
+
+        df = df.rename_geometry("geom2")
+        assert df.crs == self.osgb
+        assert df.geometry.crs == self.osgb
+        assert df.geometry.values.crs == self.osgb
+
+        df = df.rename(columns={"col1": "column1"})
+        assert df.column1.crs == self.wgs
+        assert df.column1.values.crs == self.wgs
+
+    def test_to_crs(self):
+        s = GeoSeries(self.geoms, crs=27700)
+        s = s.to_crs(4326)
+        assert s.crs == self.wgs
+        assert s.values.crs == self.wgs
+
+        df = GeoDataFrame(geometry=s)
+        assert df.crs == self.wgs
+        df = df.to_crs(27700)
+        assert df.crs == self.osgb
+        assert df.geometry.crs == self.osgb
+        assert df.geometry.values.crs == self.osgb
+
+        # make sure that only active geometry is transformed
+        arr = from_shapely(self.geoms, crs=4326)
+        df["col1"] = arr
+        df = df.to_crs(3857)
+        assert df.col1.crs == self.wgs
+        assert df.col1.values.crs == self.wgs
+
+    def test_from_shapely(self):
+        arr = from_shapely(self.geoms, crs=27700)
+        assert arr.crs == self.osgb
+
+    def test_from_wkb(self):
+        L_wkb = [p.wkb for p in self.geoms]
+        arr = from_wkb(L_wkb, crs=27700)
+        assert arr.crs == self.osgb
+
+    def test_from_wkt(self):
+        L_wkt = [p.wkt for p in self.geoms]
+        arr = from_wkt(L_wkt, crs=27700)
+        assert arr.crs == self.osgb
+
+    def test_points_from_xy(self):
+        df = pd.DataFrame([{"x": x, "y": x, "z": x} for x in range(10)])
+        arr = points_from_xy(df["x"], df["y"], crs=27700)
+        assert arr.crs == self.osgb
+
+    # setting CRS in GeoSeries should not set it in passed array without CRS
+    def test_original(self):
+        arr = from_shapely(self.geoms)
+        s = GeoSeries(arr, crs=27700)
+        assert arr.crs is None
+        assert s.crs == self.osgb
+
+    def test_ops(self):
+        arr = self.arr
+        bound = arr.boundary
+        assert bound.crs == self.osgb
+
+        cent = arr.centroid
+        assert cent.crs == self.osgb
+
+        hull = arr.convex_hull
+        assert hull.crs == self.osgb
+
+        envelope = arr.envelope
+        assert envelope.crs == self.osgb
+
+        exterior = arr.exterior
+        assert exterior.crs == self.osgb
+
+        representative_point = arr.representative_point()
+        assert representative_point.crs == self.osgb
+
+    def test_binary_ops(self):
+        arr = self.arr
+        quads = []
+        while len(quads) < 10:
+            geom = Polygon([(random.random(), random.random()) for i in range(4)])
+            if geom.is_valid:
+                quads.append(geom)
+
+        arr2 = from_shapely(quads, crs=27700)
+
+        difference = arr.difference(arr2)
+        assert difference.crs == self.osgb
+
+        intersection = arr.intersection(arr2)
+        assert intersection.crs == self.osgb
+
+        symmetric_difference = arr.symmetric_difference(arr2)
+        assert symmetric_difference.crs == self.osgb
+
+        union = arr.union(arr2)
+        assert union.crs == self.osgb
+
+    def test_other(self):
+        arr = self.arr
+
+        buffer = arr.buffer(5)
+        assert buffer.crs == self.osgb
+
+        interpolate = arr.exterior.interpolate(0.1)
+        assert interpolate.crs == self.osgb
+
+        simplify = arr.simplify(5)
+        assert simplify.crs == self.osgb
+
+    @pytest.mark.parametrize(
+        "attr, arg",
+        [
+            ("affine_transform", ([0, 1, 1, 0, 0, 0],)),
+            ("translate", ()),
+            ("rotate", (10,)),
+            ("scale", ()),
+            ("skew", ()),
+        ],
+    )
+    def test_affinity_methods(self, attr, arg):
+        result = getattr(self.arr, attr)(*arg)
+
+        assert result.crs == self.osgb
+
+    def test_slice(self):
+        s = GeoSeries(self.arr, crs=27700)
+        assert s.iloc[1:].values.crs == self.osgb
+
+        df = GeoDataFrame({"col1": self.arr}, geometry=s)
+        assert df.iloc[1:].geometry.values.crs == self.osgb
+        assert df.iloc[1:].col1.values.crs == self.osgb
+
+    def test_concat(self):
+        s = GeoSeries(self.arr, crs=27700)
+        assert pd.concat([s, s]).values.crs == self.osgb
+
+        df = GeoDataFrame({"col1": from_shapely(self.geoms, crs=4326)}, geometry=s)
+        assert pd.concat([df, df]).geometry.values.crs == self.osgb
+        assert pd.concat([df, df]).col1.values.crs == self.wgs
+
+    def test_merge(self):
+        arr = from_shapely(self.geoms, crs=27700)
+        s = GeoSeries(self.geoms, crs=4326)
+        df = GeoDataFrame({"col1": s}, geometry=arr)
+        df2 = GeoDataFrame({"col2": s}, geometry=arr).rename_geometry("geom")
+        merged = df.merge(df2, left_index=True, right_index=True)
+        assert merged.col1.values.crs == self.wgs
+        assert merged.geometry.values.crs == self.osgb
+        assert merged.col2.values.crs == self.wgs
+        assert merged.geom.values.crs == self.osgb
+        assert merged.crs == self.osgb
+
+    # CRS should be assigned to geometry
+    def test_deprecation(self):
+        with pytest.warns(FutureWarning):
+            GeoDataFrame([], crs=27700)
+
+        with pytest.warns(FutureWarning):
+            df = GeoDataFrame([])
+            df.crs = 27700
+
+    # make sure that geometry column from list has CRS (__setitem__)
+    def test_setitem_geometry(self):
+        arr = from_shapely(self.geoms, crs=27700)
+        df = GeoDataFrame({"col1": [0, 1]}, geometry=arr)
+
+        df["geometry"] = [g for g in df.geometry]
+        assert df.geometry.values.crs == self.osgb
+
+        df2 = GeoDataFrame({"col1": [0, 1]}, geometry=arr)
+        df2["geometry"] = from_shapely(self.geoms, crs=4326)
+        assert df2.geometry.values.crs == self.wgs
+
+    def test_astype(self):
+        arr = from_shapely(self.geoms, crs=27700)
+        df = GeoDataFrame({"col1": [0, 1]}, geometry=arr)
+        df2 = df.astype({"col1": str})
+        assert df2.crs == self.osgb
diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py
index 9645e7e2..c89174dc 100644
--- a/geopandas/tests/test_geodataframe.py
+++ b/geopandas/tests/test_geodataframe.py
@@ -165,7 +165,7 @@ class TestDataFrame:
         assert_geoseries_equal(df["geometry"], new_geom)
 
         # new crs
-        gs = GeoSeries(new_geom, crs="epsg:3857")
+        gs = new_geom.to_crs(crs="epsg:3857")
         df.geometry = gs
         assert df.crs == "epsg:3857"
 
@@ -214,7 +214,7 @@ class TestDataFrame:
 
         df2 = self.df.set_geometry(geom)
         assert self.df is not df2
-        assert_geoseries_equal(df2.geometry, geom)
+        assert_geoseries_equal(df2.geometry, geom, check_crs=False)
         assert_geoseries_equal(self.df.geometry, original_geom)
         assert_geoseries_equal(self.df["geometry"], self.df.geometry)
         # unknown column
diff --git a/geopandas/tests/test_testing.py b/geopandas/tests/test_testing.py
index 143faacc..b5a483e7 100644
--- a/geopandas/tests/test_testing.py
+++ b/geopandas/tests/test_testing.py
@@ -42,6 +42,19 @@ s4 = Series(a)
 df1 = GeoDataFrame({"col1": [1, 2], "geometry": s1})
 df2 = GeoDataFrame({"col1": [1, 2], "geometry": s2})
 
+s4 = s1.copy()
+s4.crs = 4326
+s5 = s2.copy()
+s5.crs = 27700
+df4 = GeoDataFrame(
+    {"col1": [1, 2], "geometry": s1.copy(), "geom2": s4.copy(), "geom3": s5.copy()},
+    crs=3857,
+)
+df5 = GeoDataFrame(
+    {"col1": [1, 2], "geometry": s1.copy(), "geom3": s5.copy(), "geom2": s4.copy()},
+    crs=3857,
+)
+
 
 def test_geoseries():
     assert_geoseries_equal(s1, s2)
@@ -69,6 +82,11 @@ def test_geodataframe():
     with pytest.raises(AssertionError):
         assert_geodataframe_equal(df1, df3)
 
+    assert_geodataframe_equal(df5, df4, check_like=True)
+    df5.geom2.crs = 3857
+    with pytest.raises(AssertionError):
+        assert_geodataframe_equal(df5, df4, check_like=True)
+
 
 def test_equal_nans():
     s = GeoSeries([Point(0, 0), np.nan])
@@ -83,7 +101,7 @@ def test_no_crs():
 
 
 def test_ignore_crs_mismatch():
-    df1 = GeoDataFrame({"col1": [1, 2], "geometry": s1}, crs="EPSG:4326")
+    df1 = GeoDataFrame({"col1": [1, 2], "geometry": s1.copy()}, crs="EPSG:4326")
     df2 = GeoDataFrame({"col1": [1, 2], "geometry": s1}, crs="EPSG:31370")
 
     with pytest.raises(AssertionError):