Repository: geopandas/geopandas

Issue:
BUG: Plotting empty geometries generates a Type Error (shapely only)
I'm not sure if this is a bug or expected behaviour. I've tested this in geopandas 0.7 and 0.8.1.

I have a process which generates some geometry but it may fail. I use that geometry for further operations so it seems prefereable to keep the failed geometries as empty geometries (e.g. Polygon()) rathern than missing (None). However when I plot the final geodataframe the plot will fail with a TypeError if the gemetries are empty but plot as expected if they are missing.

The plotting issue is easily reproducible with the [Creating a GeoDataFrame from a DataFrame with coordinates example](https://geopandas.org/gallery/create_geopandas_from_pandas.html). I would expect it to plot the same as the version with None.

```
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point

df = pd.DataFrame(
    {'City': ['Buenos Aires', 'Brasilia', 'Santiago', 'Bogota', 'Caracas'],
     'Country': ['Argentina', 'Brazil', 'Chile', 'Colombia', 'Venezuela'],
     'Latitude': [-34.58, -15.78, -33.45, 4.60, 10.48],
     'Longitude': [-58.66, -47.91, -70.66, -74.08, -66.86]})

gdf = gpd.GeoDataFrame(
    df, geometry=gpd.points_from_xy(df.Longitude, df.Latitude))

gdf.plot(color='red') # plots okay
```
Missing geometry plots okay
```
gdf.loc[3,'geometry'] = None
gdf.plot(color='red') # plots okay
```
Empty geometry raises a TypeError
```
gdf.loc[3,'geometry'] = Point()
gdf.plot(color='red') # raises a TypeError
```

Discussion:

Patch:
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index 44a23fe1..2b4168ed 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -27,6 +27,7 @@ General methods and attributes
    GeoSeries.interiors
    GeoSeries.x
    GeoSeries.y
+   GeoSeries.z
 
 Unary predicates
 ----------------
diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py
index 9868d255..b96ddacf 100644
--- a/geopandas/_vectorized.py
+++ b/geopandas/_vectorized.py
@@ -896,6 +896,14 @@ def get_y(data):
         return _unary_op("y", data, null_value=np.nan)
 
 
+def get_z(data):
+    if compat.USE_PYGEOS:
+        return pygeos.get_z(data)
+    else:
+        data = [geom.z if geom.has_z else np.nan for geom in data]
+        return np.array(data, dtype=np.dtype(float))
+
+
 def bounds(data):
     if compat.USE_PYGEOS:
         return pygeos.bounds(data)
diff --git a/geopandas/array.py b/geopandas/array.py
index ea0eb72a..bc1e302e 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -730,6 +730,15 @@ class GeometryArray(ExtensionArray):
             message = "y attribute access only provided for Point geometries"
             raise ValueError(message)
 
+    @property
+    def z(self):
+        """Return the z location of point geometries in a GeoSeries"""
+        if (self.geom_type[~self.isna()] == "Point").all():
+            return vectorized.get_z(self.data)
+        else:
+            message = "z attribute access only provided for Point geometries"
+            raise ValueError(message)
+
     @property
     def bounds(self):
         return vectorized.bounds(self.data)
diff --git a/geopandas/geoseries.py b/geopandas/geoseries.py
index 793cf7c9..82c6bd92 100644
--- a/geopandas/geoseries.py
+++ b/geopandas/geoseries.py
@@ -257,6 +257,7 @@ class GeoSeries(GeoPandasBase, Series):
         --------
 
         GeoSeries.y
+        GeoSeries.z
 
         """
         return _delegate_property("x", self)
@@ -284,10 +285,39 @@ class GeoSeries(GeoPandasBase, Series):
         --------
 
         GeoSeries.x
+        GeoSeries.z
 
         """
         return _delegate_property("y", self)
 
+    @property
+    def z(self):
+        """Return the z location of point geometries in a GeoSeries
+
+        Returns
+        -------
+        pandas.Series
+
+        Examples
+        --------
+
+        >>> from shapely.geometry import Point
+        >>> s = geopandas.GeoSeries([Point(1, 1, 1), Point(2, 2, 2), Point(3, 3, 3)])
+        >>> s.z
+        0    1.0
+        1    2.0
+        2    3.0
+        dtype: float64
+
+        See Also
+        --------
+
+        GeoSeries.x
+        GeoSeries.y
+
+        """
+        return _delegate_property("z", self)
+
     @classmethod
     def from_file(cls, filename, **kwargs):
         """Alternate constructor to create a ``GeoSeries`` from a file.
diff --git a/geopandas/plotting.py b/geopandas/plotting.py
index 275175e7..a21e9487 100644
--- a/geopandas/plotting.py
+++ b/geopandas/plotting.py
@@ -45,7 +45,7 @@ def _flatten_multi_geoms(geoms, prefix="Multi"):
         return geoms, np.arange(len(geoms))
 
     for ix, geom in enumerate(geoms):
-        if geom.type.startswith(prefix) and not geom.is_empty:
+        if geom is not None and geom.type.startswith(prefix) and not geom.is_empty:
             for poly in geom.geoms:
                 components.append(poly)
                 component_index.append(ix)
@@ -283,8 +283,8 @@ def _plot_point_collection(
     geoms, multiindex = _flatten_multi_geoms(geoms)
     # values are expanded below as kwargs["c"]
 
-    x = [p.x for p in geoms]
-    y = [p.y for p in geoms]
+    x = [p.x if not p.is_empty else None for p in geoms]
+    y = [p.y if not p.is_empty else None for p in geoms]
 
     # matplotlib 1.4 does not support c=None, and < 2.0 does not support s=None
     if values is not None:


Test patch:
diff --git a/geopandas/tests/test_array.py b/geopandas/tests/test_array.py
index 9b0f0e0d..4cc6e687 100644
--- a/geopandas/tests/test_array.py
+++ b/geopandas/tests/test_array.py
@@ -117,9 +117,11 @@ def test_from_shapely_geo_interface():
             return {"type": "Point", "coordinates": (self.x, self.y)}
 
     result = from_shapely([Point(1.0, 2.0), Point(3.0, 4.0)])
+
     expected = from_shapely(
         [shapely.geometry.Point(1.0, 2.0), shapely.geometry.Point(3.0, 4.0)]
     )
+
     assert all(v.equals(t) for v, t in zip(result, expected))
 
 
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index e4940424..1764a2f7 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -30,7 +30,10 @@ class TestGeomMethods:
         self.t1 = Polygon([(0, 0), (1, 0), (1, 1)])
         self.t2 = Polygon([(0, 0), (1, 1), (0, 1)])
         self.t3 = Polygon([(2, 0), (3, 0), (3, 1)])
+        self.tz = Polygon([(1, 1, 1), (2, 2, 2), (3, 3, 3)])
+        self.tz1 = Polygon([(2, 2, 2), (1, 1, 1), (3, 3, 3)])
         self.sq = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
+        self.sqz = Polygon([(1, 1, 1), (2, 2, 2), (3, 3, 3), (4, 4, 4)])
         self.t4 = Polygon([(0, 0), (3, 0), (3, 3), (0, 2)])
         self.t5 = Polygon([(2, 0), (3, 0), (3, 3), (2, 3)])
         self.inner_sq = Polygon(
@@ -53,6 +56,7 @@ class TestGeomMethods:
         self.g1 = GeoSeries([self.t1, self.sq])
         self.g2 = GeoSeries([self.sq, self.t1])
         self.g3 = GeoSeries([self.t1, self.t2])
+        self.gz = GeoSeries([self.tz, self.sqz, self.tz1])
         self.g3.crs = "epsg:4326"
         self.g4 = GeoSeries([self.t2, self.t1])
         self.g4.crs = "epsg:4326"
@@ -63,9 +67,11 @@ class TestGeomMethods:
         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.esb = Point(-73.9847, 40.7484, 30.3244)
+        self.sol = Point(-74.0446, 40.6893, 31.2344)
         self.landmarks = GeoSeries([self.esb, self.sol], crs="epsg:4326")
+        self.pt2d = Point(-73.9847, 40.7484)
+        self.landmarks_mixed = GeoSeries([self.esb, self.sol, self.pt2d], crs=4326)
         self.l1 = LineString([(0, 0), (0, 1), (1, 1)])
         self.l2 = LineString([(0, 0), (1, 0), (1, 1), (0, 1)])
         self.g5 = GeoSeries([self.l1, self.l2])
@@ -92,6 +98,9 @@ class TestGeomMethods:
         self.gdf3 = GeoDataFrame(
             {"geometry": self.g3, "col3": [4, 5], "col4": ["rand", "string"]}
         )
+        self.gdfz = GeoDataFrame(
+            {"geometry": self.gz, "col3": [4, 5, 6], "col4": ["rand", "string", "geo"]}
+        )
 
     def _test_unary_real(self, op, expected, a):
         """ Tests for 'area', 'length', 'is_valid', etc. """
@@ -464,20 +473,29 @@ class TestGeomMethods:
         expected = Series([False, True], self.g_3d.index)
         self._test_unary_real("has_z", expected, self.g_3d)
 
-    def test_xy_points(self):
+    def test_xyz_points(self):
         expected_x = [-73.9847, -74.0446]
         expected_y = [40.7484, 40.6893]
+        expected_z = [30.3244, 31.2344]
 
         assert_array_dtype_equal(expected_x, self.landmarks.geometry.x)
         assert_array_dtype_equal(expected_y, self.landmarks.geometry.y)
+        assert_array_dtype_equal(expected_z, self.landmarks.geometry.z)
+
+        # mixed dimensions
+        expected_z = [30.3244, 31.2344, np.nan]
+        assert_array_dtype_equal(expected_z, self.landmarks_mixed.geometry.z)
 
-    def test_xy_polygons(self):
+    def test_xyz_polygons(self):
         # accessing x attribute in polygon geoseries should raise an error
         with pytest.raises(ValueError):
             _ = self.gdf1.geometry.x
         # and same for accessing y attribute in polygon geoseries
         with pytest.raises(ValueError):
             _ = self.gdf1.geometry.y
+        # and same for accessing z attribute in polygon geoseries
+        with pytest.raises(ValueError):
+            _ = self.gdfz.geometry.z
 
     def test_centroid(self):
         polygon = Polygon([(-1, -1), (1, -1), (1, 1), (-1, 1)])
diff --git a/geopandas/tests/test_plotting.py b/geopandas/tests/test_plotting.py
index cb518550..ccee1083 100644
--- a/geopandas/tests/test_plotting.py
+++ b/geopandas/tests/test_plotting.py
@@ -308,6 +308,18 @@ class TestPointPlotting:
             ax = s.plot()
             assert len(ax.collections) == 1
 
+        # more complex case with GEOMETRYCOLLECTION EMPTY, POINT EMPTY and NONE
+        poly = Polygon([(-1, -1), (-1, 2), (2, 2), (2, -1), (-1, -1)])
+        point = Point(0, 1)
+        point_ = Point(10, 10)
+        empty_point = Point()
+
+        gdf = GeoDataFrame(geometry=[point, empty_point, point_])
+        gdf["geometry"] = gdf.intersection(poly)
+        gdf.loc[3] = [None]
+        ax = gdf.plot()
+        assert len(ax.collections) == 1
+
     def test_multipoints(self):
 
         # MultiPoints