Repository: geopandas/geopandas
Pull request #3133 (2024-01-06 14:56:05)

ENH: harmonize read_file column selection keywords (support `columns` for fiona, support `ignore_fields`/`include_fields` for pyogrio for back compat)
A couple of changes related to selecting columns in `read_file`:

- Publicize to use `columns` (which I think the nicest API): explicitly add to the docstring (and will update the docs to mention it as well)
- Translate `columns` to `incude_fields` for the fiona engine, so `columns` is always supported.
- For the pyogrio engine, translate both `include_fields` and `ignore_fields` keywords that are supported by fiona to the `columns` keyword of pyogrio. This ensures a smooth upgrade when we switch the default for people that were using those keywords.
  - I did add a warning pointing to use `columns` instead, we can discuss whether we want to keep this warning (or make it stronger).



martinfleis: +1 on keeping only `columns` in future. 
jorisvandenbossche: Updated the PR. Fixed the failing test (remaining failure is the dev build, failing for other reasons). And changed the UserWarning into a DeprecationWarning with the intent that eventually we only keep the `columns` keyword.

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7393cf5a..d4539b99 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,6 +11,8 @@ Notes on dependencies:
 API changes:
 
 - `unary_union` is now deprecated and replaced by the `union_all` method (#3007).
+- `align` keyword in binary methods now defaults to `None`, treated as True. Explicit True
+  will silence the warning about mismachted indices. (#3212)
 - The `sjoin` method will now preserve the name of the index of the right
   GeoDataFrame, if it has one, instead of always using `"index_right"` as the
   name for the resulting column in the return value (#846, #2144).
@@ -42,9 +44,8 @@ New features and improvements:
   `buffer(0)` (#3113).
 - Passing `"geometry"` as `dtype` to `pd.read_csv` will now return a GeoSeries for
   the specified columns (#3101).
--
-API changes:
-- Added support for ``mask`` keyword for pyogrio engine for pyogrio >= 0.7.0 (#3062).
+- Added support to ``read_file`` for the ``mask`` keyword for the pyogrio engine (#3062).
+- Added support to ``read_file`` for the ``columns`` keyword for the fiona engine (#3133).
 
 Backwards incompatible API changes:
 - The deprecated default value of GeoDataFrame/ GeoSeries `explode(.., index_parts=True)` is now
@@ -80,7 +81,10 @@ Deprecations and compatibility notes:
     `__sub__`. Instead use methods `symmetric_difference`, `union`, `intersection` and
     `difference` respectively.
 - Fixes for compatibility with psycopg (#3167).
-
+- The ``include_fields`` and ``ignore_fields`` keywords in ``read_file()`` are deprecated
+  for the default pyogrio engine. Currently those are translated to the ``columns`` keyword
+  for backwards compatibility, but you should directly use the ``columns`` keyword instead
+  to select which columns to read (#3133).
 
 ## Version 0.14.3 (Jan 31, 2024)
 
diff --git a/ci/envs/311-dev.yaml b/ci/envs/311-dev.yaml
index bc7964a3..313e8758 100644
--- a/ci/envs/311-dev.yaml
+++ b/ci/envs/311-dev.yaml
@@ -5,7 +5,6 @@ dependencies:
   - python=3.11
   - cython
   # required
-  - numpy
   - pyproj
   - geos
   - packaging
@@ -24,6 +23,7 @@ dependencies:
     - mapclassify>=2.4.0
     # dev versions of packages
     - --pre --index-url https://pypi.anaconda.org/scientific-python-nightly-wheels/simple --extra-index-url https://pypi.fury.io/arrow-nightlies/ --extra-index-url https://pypi.org/simple
+    - numpy>=2.0.0.dev
     - fiona
     - pandas
     - matplotlib
diff --git a/doc/source/docs/user_guide/io.rst b/doc/source/docs/user_guide/io.rst
index 3c597898..58aadb10 100644
--- a/doc/source/docs/user_guide/io.rst
+++ b/doc/source/docs/user_guide/io.rst
@@ -143,29 +143,18 @@ or a slice object.
 Field/column filters
 ^^^^^^^^^^^^^^^^^^^^
 
-Load in a subset of fields from the file:
-
-.. note:: Requires Fiona 1.9+
-
-.. code-block:: python
-
-    gdf = geopandas.read_file(
-        geodatasets.get_path("geoda.nyc"),
-        include_fields=["name", "rent2008", "kids2000"],
-    )
-
-.. note:: Requires Fiona 1.8+
+Load in a subset of fields from the file using the ``columns`` keyword
+(this requires pyogrio or Fiona 1.9+):
 
 .. code-block:: python
 
     gdf = geopandas.read_file(
         geodatasets.get_path("geoda.nyc"),
-        ignore_fields=["rent2008", "kids2000"],
+        columns=["name", "rent2008", "kids2000"],
     )
 
 Skip loading geometry from the file:
 
-.. note:: Requires Fiona 1.8+
 .. note:: Returns :obj:`pandas.DataFrame`
 
 .. code-block:: python
diff --git a/geopandas/base.py b/geopandas/base.py
index 4691fe1e..c19b4d25 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -27,13 +27,24 @@ def is_geometry_type(data):
 
 def _delegate_binary_method(op, this, other, align, *args, **kwargs):
     # type: (str, GeoSeries, GeoSeries) -> GeoSeries/Series
+    if align is None:
+        align = True
+        maybe_warn = True
+    else:
+        maybe_warn = False
     this = this.geometry
     if isinstance(other, GeoPandasBase):
         if align and not this.index.equals(other.index):
-            warn(
-                "The indices of the two GeoSeries are different.",
-                stacklevel=4,
-            )
+            if maybe_warn:
+                warn(
+                    "The indices of the left and right GeoSeries' are not equal, and "
+                    "therefore they will be aligned (reordering and/or introducing "
+                    "missing values) before executing the operation. If this alignment "
+                    "is the desired behaviour, you can silence this warning by passing "
+                    "'align=True'. If you don't want alignment and protect yourself of "
+                    "accidentally aligning, you can pass 'align=False'.",
+                    stacklevel=4,
+                )
             this, other = this.align(other.geometry)
         else:
             other = other.geometry
@@ -1737,7 +1748,7 @@ GeometryCollection
     # Binary operations that return a pandas Series
     #
 
-    def contains(self, other, align=True):
+    def contains(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that contains `other`.
 
@@ -1760,9 +1771,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if it
             is contained.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -1853,7 +1864,7 @@ GeometryCollection
         """
         return _binary_op("contains", self, other, align)
 
-    def contains_properly(self, other, align=True):
+    def contains_properly(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that is completely inside ``other``, with no common
         boundary points.
@@ -1873,9 +1884,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if it
             is contained.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -1974,7 +1985,7 @@ GeometryCollection
         """
         return _binary_op("contains_properly", self, other, align)
 
-    def dwithin(self, other, distance, align=True):
+    def dwithin(self, other, distance, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that is within a set distance from ``other``.
 
@@ -1993,9 +2004,9 @@ GeometryCollection
             applied to all geometries. An array or Series will be applied elementwise.
             If np.array or pd.Series are used then it must have same length as the
             GeoSeries.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2085,7 +2096,7 @@ GeometryCollection
         """
         return _binary_op("dwithin", self, other, distance=distance, align=align)
 
-    def geom_equals(self, other, align=True):
+    def geom_equals(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry equal to `other`.
 
@@ -2103,9 +2114,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test for
             equality.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2195,7 +2206,7 @@ GeometryCollection
         """
         return _binary_op("geom_equals", self, other, align)
 
-    def geom_almost_equals(self, other, decimal=6, align=True):
+    def geom_almost_equals(self, other, decimal=6, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` if
         each aligned geometry is approximately equal to `other`.
 
@@ -2213,9 +2224,9 @@ GeometryCollection
             The GeoSeries (elementwise) or geometric object to compare to.
         decimal : int
             Decimal place precision used when testing for approximate equality.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2272,7 +2283,7 @@ GeometryCollection
             "geom_equals_exact", self, other, tolerance=tolerance, align=align
         )
 
-    def geom_equals_exact(self, other, tolerance, align=True):
+    def geom_equals_exact(self, other, tolerance, align=None):
         """Return True for all geometries that equal aligned *other* to a given
         tolerance, else False.
 
@@ -2287,9 +2298,9 @@ GeometryCollection
             The GeoSeries (elementwise) or geometric object to compare to.
         tolerance : float
             Decimal place precision used when testing for approximate equality.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2337,7 +2348,7 @@ GeometryCollection
             "geom_equals_exact", self, other, tolerance=tolerance, align=align
         )
 
-    def crosses(self, other, align=True):
+    def crosses(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that cross `other`.
 
@@ -2355,9 +2366,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if is
             crossed.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2449,7 +2460,7 @@ GeometryCollection
         """
         return _binary_op("crosses", self, other, align)
 
-    def disjoint(self, other, align=True):
+    def disjoint(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry disjoint to `other`.
 
@@ -2466,9 +2477,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if is
             disjoint.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2550,7 +2561,7 @@ GeometryCollection
         """
         return _binary_op("disjoint", self, other, align)
 
-    def intersects(self, other, align=True):
+    def intersects(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that intersects `other`.
 
@@ -2567,9 +2578,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if is
             intersected.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2661,7 +2672,7 @@ GeometryCollection
         """
         return _binary_op("intersects", self, other, align)
 
-    def overlaps(self, other, align=True):
+    def overlaps(self, other, align=None):
         """Returns True for all aligned geometries that overlap *other*, else False.
 
         Geometries overlaps if they have more than one but not all
@@ -2679,9 +2690,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if
             overlaps.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2772,7 +2783,7 @@ GeometryCollection
         """
         return _binary_op("overlaps", self, other, align)
 
-    def touches(self, other, align=True):
+    def touches(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that touches `other`.
 
@@ -2790,9 +2801,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if is
             touched.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2884,7 +2895,7 @@ GeometryCollection
         """
         return _binary_op("touches", self, other, align)
 
-    def within(self, other, align=True):
+    def within(self, other, align=None):
         """Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that is within `other`.
 
@@ -2906,9 +2917,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The GeoSeries (elementwise) or geometric object to test if each
             geometry is within.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -2998,7 +3009,7 @@ GeometryCollection
         """
         return _binary_op("within", self, other, align)
 
-    def covers(self, other, align=True):
+    def covers(self, other, align=None):
         """
         Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that is entirely covering `other`.
@@ -3020,9 +3031,9 @@ GeometryCollection
         ----------
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to check is being covered.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -3112,7 +3123,7 @@ GeometryCollection
         """
         return _binary_op("covers", self, other, align)
 
-    def covered_by(self, other, align=True):
+    def covered_by(self, other, align=None):
         """
         Returns a ``Series`` of ``dtype('bool')`` with value ``True`` for
         each aligned geometry that is entirely covered by `other`.
@@ -3133,9 +3144,9 @@ GeometryCollection
         ----------
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to check is being covered.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -3226,7 +3237,7 @@ GeometryCollection
         """
         return _binary_op("covered_by", self, other, align)
 
-    def distance(self, other, align=True):
+    def distance(self, other, align=None):
         """Returns a ``Series`` containing the distance to aligned `other`.
 
         The operation works on a 1-to-1 row-wise manner:
@@ -3239,9 +3250,9 @@ GeometryCollection
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             distance to.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
 
         Returns
@@ -3338,9 +3349,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             distance to.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
         densify : float (default None)
             A value between 0 and 1, that splits each subsegment of a line string
             into equal length segments, making the approximation less coarse.
@@ -3460,9 +3471,9 @@ GeometryCollection
         other : GeoSeries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             distance to.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
         densify : float (default None)
             A value between 0 and 1, that splits each subsegment of a line string
             into equal length segments, making the approximation less coarse.
@@ -3561,7 +3572,7 @@ GeometryCollection
     # Binary operations that return a GeoSeries
     #
 
-    def difference(self, other, align=True):
+    def difference(self, other, align=None):
         """Returns a ``GeoSeries`` of the points in each aligned geometry that
         are not in `other`.
 
@@ -3578,9 +3589,9 @@ GeometryCollection
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             difference to.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -3672,7 +3683,7 @@ GeometryCollection
         """
         return _binary_geo("difference", self, other, align)
 
-    def symmetric_difference(self, other, align=True):
+    def symmetric_difference(self, other, align=None):
         """Returns a ``GeoSeries`` of the symmetric difference of points in
         each aligned geometry with `other`.
 
@@ -3693,9 +3704,9 @@ GeometryCollection
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             symmetric difference to.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -3787,7 +3798,7 @@ GeometryCollection
         """
         return _binary_geo("symmetric_difference", self, other, align)
 
-    def union(self, other, align=True):
+    def union(self, other, align=None):
         """Returns a ``GeoSeries`` of the union of points in each aligned geometry with
         `other`.
 
@@ -3805,9 +3816,9 @@ GeometryCollection
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to find the union
             with.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -3901,7 +3912,7 @@ GeometryCollection
         """
         return _binary_geo("union", self, other, align)
 
-    def intersection(self, other, align=True):
+    def intersection(self, other, align=None):
         """Returns a ``GeoSeries`` of the intersection of points in each
         aligned geometry with `other`.
 
@@ -3919,9 +3930,9 @@ GeometryCollection
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             intersection with.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -4083,7 +4094,7 @@ GeometryCollection
         clipped_geometry = geometry_array.clip_by_rect(xmin, ymin, xmax, ymax)
         return GeoSeries(clipped_geometry, index=self.index, crs=self.crs)
 
-    def shortest_line(self, other, align=True):
+    def shortest_line(self, other, align=None):
         """
         Returns the shortest two-point line between two geometries.
 
@@ -4104,9 +4115,9 @@ GeometryCollection
         other : Geoseries or geometric object
             The Geoseries (elementwise) or geometric object to find the
             shortest line with.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -4186,7 +4197,7 @@ GeometryCollection
         """
         return _binary_geo("shortest_line", self, other, align)
 
-    def snap(self, other, tolerance, align=True):
+    def snap(self, other, tolerance, align=None):
         """Snaps an input geometry to reference geometry’s vertices.
 
         Vertices of the first geometry are snapped to vertices of the second. geometry,
@@ -4215,9 +4226,9 @@ GeometryCollection
             The Geoseries (elementwise) or geometric object to snap to.
         tolerance : float or array like
             Maximum distance between vertices that shall be snapped
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -4572,7 +4583,7 @@ GeometryCollection
         """
         return _delegate_geo_method("simplify", self, *args, **kwargs)
 
-    def relate(self, other, align=True):
+    def relate(self, other, align=None):
         """
         Returns the DE-9IM intersection matrices for the geometries
 
@@ -4586,9 +4597,9 @@ GeometryCollection
         other : BaseGeometry or GeoSeries
             The other geometry to computed
             the DE-9IM intersection matrices from.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
@@ -4677,7 +4688,7 @@ GeometryCollection
         """
         return _binary_op("relate", self, other, align)
 
-    def project(self, other, normalized=False, align=True):
+    def project(self, other, normalized=False, align=None):
         """
         Return the distance along each geometry nearest to *other*
 
@@ -4698,9 +4709,9 @@ GeometryCollection
         normalized : boolean
             If normalized is True, return the distance normalized to
             the length of the object.
-        align : bool (default True)
+        align : bool | None (default None)
             If True, automatically aligns GeoSeries based on their indices.
-            If False, the order of elements is preserved.
+            If False, the order of elements is preserved. None defaults to True.
 
         Returns
         -------
diff --git a/geopandas/io/file.py b/geopandas/io/file.py
index a3cd6ddc..74bcd894 100644
--- a/geopandas/io/file.py
+++ b/geopandas/io/file.py
@@ -177,7 +177,9 @@ def _is_url(url):
         return False
 
 
-def _read_file(filename, bbox=None, mask=None, rows=None, engine=None, **kwargs):
+def _read_file(
+    filename, bbox=None, mask=None, columns=None, rows=None, engine=None, **kwargs
+):
     """
     Returns a GeoDataFrame from a file or URL.
 
@@ -210,6 +212,11 @@ def _read_file(filename, bbox=None, mask=None, rows=None, engine=None, **kwargs)
         CRS mis-matches are resolved if given a GeoSeries or GeoDataFrame.
         Cannot be used with bbox. If multiple geometries are passed, this will
         first union all geometries, which may be computationally expensive.
+    columns : list, optional
+        List of column names to import from the data source. Column names
+        must exactly match the names in the data source. To avoid reading
+        any columns (besides the geometry column), pass an empty list-like.
+        By default reads all columns.
     rows : int or slice, default None
         Load in specific rows by passing an integer (first `n` rows) or a
         slice() object.
@@ -284,7 +291,9 @@ def _read_file(filename, bbox=None, mask=None, rows=None, engine=None, **kwargs)
                 from_bytes = True
 
     if engine == "pyogrio":
-        return _read_file_pyogrio(filename, bbox=bbox, mask=mask, rows=rows, **kwargs)
+        return _read_file_pyogrio(
+            filename, bbox=bbox, mask=mask, columns=columns, rows=rows, **kwargs
+        )
 
     elif engine == "fiona":
         if pd.api.types.is_file_like(filename):
@@ -295,7 +304,13 @@ def _read_file(filename, bbox=None, mask=None, rows=None, engine=None, **kwargs)
             path_or_bytes = filename
 
         return _read_file_fiona(
-            path_or_bytes, from_bytes, bbox=bbox, mask=mask, rows=rows, **kwargs
+            path_or_bytes,
+            from_bytes,
+            bbox=bbox,
+            mask=mask,
+            columns=columns,
+            rows=rows,
+            **kwargs,
         )
 
     else:
@@ -303,11 +318,31 @@ def _read_file(filename, bbox=None, mask=None, rows=None, engine=None, **kwargs)
 
 
 def _read_file_fiona(
-    path_or_bytes, from_bytes, bbox=None, mask=None, rows=None, where=None, **kwargs
+    path_or_bytes,
+    from_bytes,
+    bbox=None,
+    mask=None,
+    columns=None,
+    rows=None,
+    where=None,
+    **kwargs,
 ):
     if where is not None and not FIONA_GE_19:
         raise NotImplementedError("where requires fiona 1.9+")
 
+    if columns is not None:
+        if "include_fields" in kwargs:
+            raise ValueError(
+                "Cannot specify both 'include_fields' and 'columns' keywords"
+            )
+        if not FIONA_GE_19:
+            raise NotImplementedError("'columns' keyword requires fiona 1.9+")
+        kwargs["include_fields"] = columns
+    elif "include_fields" in kwargs:
+        # alias to columns, as this variable is used below to specify column order
+        # in the dataframe creation
+        columns = kwargs["include_fields"]
+
     if not from_bytes:
         # Opening a file via URL or file-like-object above automatically detects a
         # zipped file. In order to match that behavior, attempt to add a zip scheme
@@ -368,7 +403,7 @@ def _read_file_fiona(
             else:
                 f_filt = features
             # get list of columns
-            columns = list(features.schema["properties"])
+            columns = columns or list(features.schema["properties"])
             datetime_fields = [
                 k for (k, v) in features.schema["properties"].items() if v == "datetime"
             ]
@@ -480,6 +515,40 @@ def _read_file_pyogrio(path_or_bytes, bbox=None, mask=None, rows=None, **kwargs)
     if kwargs.pop("ignore_geometry", False):
         kwargs["read_geometry"] = False
 
+    # translate `ignore_fields`/`include_fields` keyword for back compat with fiona
+    if "ignore_fields" in kwargs and "include_fields" in kwargs:
+        raise ValueError("Cannot specify both 'ignore_fields' and 'include_fields'")
+    elif "ignore_fields" in kwargs:
+        if kwargs.get("columns", None) is not None:
+            raise ValueError(
+                "Cannot specify both 'columns' and 'ignore_fields' keywords"
+            )
+        warnings.warn(
+            "The 'include_fields' and 'ignore_fields' keywords are deprecated, and "
+            "will be removed in a future release. You can use the 'columns' keyword "
+            "instead to select which columns to read.",
+            DeprecationWarning,
+            stacklevel=3,
+        )
+        ignore_fields = kwargs.pop("ignore_fields")
+        fields = pyogrio.read_info(path_or_bytes)["fields"]
+        include_fields = [col for col in fields if col not in ignore_fields]
+        kwargs["columns"] = include_fields
+    elif "include_fields" in kwargs:
+        # translate `include_fields` keyword for back compat with fiona engine
+        if kwargs.get("columns", None) is not None:
+            raise ValueError(
+                "Cannot specify both 'columns' and 'include_fields' keywords"
+            )
+        warnings.warn(
+            "The 'include_fields' and 'ignore_fields' keywords are deprecated, and "
+            "will be removed in a future release. You can use the 'columns' keyword "
+            "instead to select which columns to read.",
+            DeprecationWarning,
+            stacklevel=3,
+        )
+        kwargs["columns"] = kwargs.pop("include_fields")
+
     return pyogrio.read_dataframe(path_or_bytes, bbox=bbox, **kwargs)
 
 


Test patch:
diff --git a/geopandas/io/tests/test_file.py b/geopandas/io/tests/test_file.py
index 351b2752..2904eb85 100644
--- a/geopandas/io/tests/test_file.py
+++ b/geopandas/io/tests/test_file.py
@@ -790,12 +790,26 @@ def test_read_file__ignore_geometry(engine, naturalearth_lowres):
     assert isinstance(pdf, pd.DataFrame) and not isinstance(pdf, geopandas.GeoDataFrame)
 
 
+@pytest.mark.filterwarnings(
+    "ignore:The 'include_fields' and 'ignore_fields' keywords:DeprecationWarning"
+)
+def test_read_file__ignore_fields(engine, naturalearth_lowres):
+    gdf = geopandas.read_file(
+        naturalearth_lowres,
+        ignore_fields=["pop_est", "continent", "iso_a3", "gdp_md_est"],
+        engine=engine,
+    )
+    assert gdf.columns.tolist() == ["name", "geometry"]
+
+
+@pytest.mark.filterwarnings(
+    "ignore:The 'include_fields' and 'ignore_fields' keywords:DeprecationWarning"
+)
 def test_read_file__ignore_all_fields(engine, naturalearth_lowres):
-    skip_pyogrio_not_supported(engine)  # pyogrio has "columns" keyword instead
     gdf = geopandas.read_file(
         naturalearth_lowres,
         ignore_fields=["pop_est", "continent", "name", "iso_a3", "gdp_md_est"],
-        engine="fiona",
+        engine=engine,
     )
     assert gdf.columns.tolist() == ["geometry"]
 
@@ -856,17 +870,65 @@ def test_read_file__where_filter(engine, naturalearth_lowres):
             )
 
 
-@PYOGRIO_MARK
-def test_read_file__columns(naturalearth_lowres):
-    # TODO: this is only support for pyogrio, but we could mimic it for fiona as well
+def test_read_file__columns(engine, naturalearth_lowres):
+    if engine == "fiona" and not FIONA_GE_19:
+        pytest.skip("columns requires fiona 1.9+")
+
     gdf = geopandas.read_file(
-        naturalearth_lowres,
-        columns=["name", "pop_est"],
-        engine="pyogrio",
+        naturalearth_lowres, columns=["name", "pop_est"], engine=engine
     )
     assert gdf.columns.tolist() == ["name", "pop_est", "geometry"]
 
 
+def test_read_file__columns_empty(engine, naturalearth_lowres):
+    if engine == "fiona" and not FIONA_GE_19:
+        pytest.skip("columns requires fiona 1.9+")
+
+    gdf = geopandas.read_file(naturalearth_lowres, columns=[], engine=engine)
+    assert gdf.columns.tolist() == ["geometry"]
+
+
+@pytest.mark.skipif(FIONA_GE_19, reason="test for fiona < 1.9")
+def test_read_file__columns_old_fiona(naturalearth_lowres):
+    with pytest.raises(NotImplementedError):
+        geopandas.read_file(
+            naturalearth_lowres, columns=["name", "pop_est"], engine="fiona"
+        )
+
+
+@pytest.mark.filterwarnings(
+    "ignore:The 'include_fields' and 'ignore_fields' keywords:DeprecationWarning"
+)
+def test_read_file__include_fields(engine, naturalearth_lowres):
+    if engine == "fiona" and not FIONA_GE_19:
+        pytest.skip("columns requires fiona 1.9+")
+
+    gdf = geopandas.read_file(
+        naturalearth_lowres, include_fields=["name", "pop_est"], engine=engine
+    )
+    assert gdf.columns.tolist() == ["name", "pop_est", "geometry"]
+
+
+@pytest.mark.skipif(not FIONA_GE_19, reason="columns requires fiona 1.9+")
+def test_read_file__columns_conflicting_keywords(engine, naturalearth_lowres):
+    path = naturalearth_lowres
+
+    with pytest.raises(ValueError, match="Cannot specify both"):
+        geopandas.read_file(
+            path, include_fields=["name"], ignore_fields=["pop_est"], engine=engine
+        )
+
+    with pytest.raises(ValueError, match="Cannot specify both"):
+        geopandas.read_file(
+            path, columns=["name"], include_fields=["pop_est"], engine=engine
+        )
+
+    with pytest.raises(ValueError, match="Cannot specify both"):
+        geopandas.read_file(
+            path, columns=["name"], ignore_fields=["pop_est"], engine=engine
+        )
+
+
 @pytest.mark.skipif(not HAS_PYPROJ, reason="pyproj not installed")
 @pytest.mark.parametrize("file_like", [False, True])
 def test_read_file_bbox_gdf(df_nybb, engine, nybb_filename, file_like):
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index 7fe0ef62..058594ef 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -310,15 +310,26 @@ class TestGeomMethods:
     #         self._test_binary_topological('intersection', self.g3,
     #                                       self.g3, no_crs_g3)
 
+    def test_alignment_warning(self):
+        with pytest.warns(
+            UserWarning,
+            match="The indices of the left and right GeoSeries' are not equal",
+        ):
+            self.g0.intersection(self.g9, align=None)
+
+        with warnings.catch_warnings(record=True) as record:
+            self.g0.intersection(self.g9, align=True)
+            self.g0.intersection(self.g9, align=False)
+
+            assert len(record) == 0
+
     def test_intersection(self):
         self._test_binary_topological("intersection", self.t1, self.g1, self.g2)
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            self._test_binary_topological(
-                "intersection", self.all_none, self.g1, self.empty
-            )
+        self._test_binary_topological(
+            "intersection", self.all_none, self.g1, self.empty
+        )
 
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert len(self.g0.intersection(self.g9, align=True) == 8)
+        assert len(self.g0.intersection(self.g9, align=True) == 8)
         assert len(self.g0.intersection(self.g9, align=False) == 7)
 
     def test_clip_by_rect(self):
@@ -333,8 +344,7 @@ class TestGeomMethods:
     def test_union_series(self):
         self._test_binary_topological("union", self.sq, self.g1, self.g2)
 
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert len(self.g0.union(self.g9, align=True) == 8)
+        assert len(self.g0.union(self.g9, align=True) == 8)
         assert len(self.g0.union(self.g9, align=False) == 7)
 
     def test_union_polygon(self):
@@ -343,8 +353,7 @@ class TestGeomMethods:
     def test_symmetric_difference_series(self):
         self._test_binary_topological("symmetric_difference", self.sq, self.g3, self.g4)
 
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert len(self.g0.symmetric_difference(self.g9, align=True) == 8)
+        assert len(self.g0.symmetric_difference(self.g9, align=True) == 8)
         assert len(self.g0.symmetric_difference(self.g9, align=False) == 7)
 
     def test_symmetric_difference_poly(self):
@@ -357,8 +366,7 @@ class TestGeomMethods:
         expected = GeoSeries([GeometryCollection(), self.t2])
         self._test_binary_topological("difference", expected, self.g1, self.g2)
 
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert len(self.g0.difference(self.g9, align=True) == 8)
+        assert len(self.g0.difference(self.g9, align=True) == 8)
         assert len(self.g0.difference(self.g9, align=False) == 7)
 
     def test_difference_poly(self):
@@ -491,8 +499,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.contains(self.t1))
 
         expected = [False, True, True, True, True, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.contains(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.contains(self.g9, align=True))
 
         expected = [False, False, True, False, False, False, False]
         assert_array_dtype_equal(expected, self.g0.contains(self.g9, align=False))
@@ -502,10 +509,9 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.contains_properly(Point(0.25, 0.25)))
 
         expected = [False, False, False, False, False, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(
-                expected, self.g0.contains_properly(self.g9, align=True)
-            )
+        assert_array_dtype_equal(
+            expected, self.g0.contains_properly(self.g9, align=True)
+        )
 
         expected = [False, False, True, False, False, False, False]
         assert_array_dtype_equal(
@@ -518,8 +524,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.dwithin(self.p0, 6))
 
         expected = [False, True, True, True, True, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.dwithin(self.g9, 1, align=True))
+        assert_array_dtype_equal(expected, self.g0.dwithin(self.g9, 1, align=True))
         expected = [True, True, True, True, False, False, False]
         assert_array_dtype_equal(expected, self.g0.dwithin(self.g9, 1, align=False))
 
@@ -550,8 +555,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.crossed_lines.crosses(self.l3))
 
         expected = [False] * 8
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.crosses(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.crosses(self.g9, align=True))
 
         expected = [False] * 7
         assert_array_dtype_equal(expected, self.g0.crosses(self.g9, align=False))
@@ -561,8 +565,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.disjoint(self.t1))
 
         expected = [False] * 8
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.disjoint(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.disjoint(self.g9, align=True))
 
         expected = [False, False, False, False, True, False, False]
         assert_array_dtype_equal(expected, self.g0.disjoint(self.g9, align=False))
@@ -599,8 +602,7 @@ class TestGeomMethods:
             index=range(8),
         )
 
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.relate(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.relate(self.g9, align=True))
 
         expected = Series(
             [
@@ -626,8 +628,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g6.distance(self.na_none))
 
         expected = Series(np.array([np.nan, 0, 0, 0, 0, 0, np.nan, np.nan]), range(8))
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.distance(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.distance(self.g9, align=True))
 
         val = self.g0.iloc[4].distance(self.g9.iloc[4])
         expected = Series(np.array([0, 0, 0, 0, val, np.nan, np.nan]), self.g0.index)
@@ -647,10 +648,9 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.na_none.hausdorff_distance(self.p0))
 
         expected = Series(np.array([np.nan, 0, 0, 0, 0, 0, np.nan, np.nan]), range(8))
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(
-                expected, self.g0.hausdorff_distance(self.g9, align=True)
-            )
+        assert_array_dtype_equal(
+            expected, self.g0.hausdorff_distance(self.g9, align=True)
+        )
 
         val_1 = self.g0.iloc[0].hausdorff_distance(self.g9.iloc[0])
         val_2 = self.g0.iloc[2].hausdorff_distance(self.g9.iloc[2])
@@ -676,10 +676,9 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.na_none.frechet_distance(self.p0))
 
         expected = Series(np.array([np.nan, 0, 0, 0, 0, 0, np.nan, np.nan]), range(8))
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(
-                expected, self.g0.frechet_distance(self.g9, align=True)
-            )
+        assert_array_dtype_equal(
+            expected, self.g0.frechet_distance(self.g9, align=True)
+        )
 
         # expected returns
         val_1 = 1.0
@@ -715,8 +714,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.intersects(self.empty_poly))
 
         expected = [False, True, True, True, True, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.intersects(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.intersects(self.g9, align=True))
 
         expected = [True, True, True, True, False, False, False]
         assert_array_dtype_equal(expected, self.g0.intersects(self.g9, align=False))
@@ -729,8 +727,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g4.overlaps(self.t1))
 
         expected = [False] * 8
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.overlaps(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.overlaps(self.g9, align=True))
 
         expected = [False] * 7
         assert_array_dtype_equal(expected, self.g0.overlaps(self.g9, align=False))
@@ -740,8 +737,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.touches(self.t1))
 
         expected = [False] * 8
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.touches(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.touches(self.g9, align=True))
 
         expected = [True, False, False, True, False, False, False]
         assert_array_dtype_equal(expected, self.g0.touches(self.g9, align=False))
@@ -754,8 +750,7 @@ class TestGeomMethods:
         assert_array_dtype_equal(expected, self.g0.within(self.sq))
 
         expected = [False, True, True, True, True, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.within(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.within(self.g9, align=True))
 
         expected = [False, True, False, False, False, False, False]
         assert_array_dtype_equal(expected, self.g0.within(self.g9, align=False))
@@ -772,8 +767,7 @@ class TestGeomMethods:
         assert_series_equal(res, exp)
 
         expected = [False, True, True, True, True, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.covers(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.covers(self.g9, align=True))
 
         expected = [False, False, True, False, False, False, False]
         assert_array_dtype_equal(expected, self.g0.covers(self.g9, align=False))
@@ -789,8 +783,7 @@ class TestGeomMethods:
         assert_series_equal(res, exp)
 
         expected = [False, True, True, True, True, True, False, False]
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_array_dtype_equal(expected, self.g0.covered_by(self.g9, align=True))
+        assert_array_dtype_equal(expected, self.g0.covered_by(self.g9, align=True))
 
         expected = [False, True, False, False, False, False, False]
         assert_array_dtype_equal(expected, self.g0.covered_by(self.g9, align=False))
@@ -1103,8 +1096,7 @@ class TestGeomMethods:
 
         s = GeoSeries([Point(2, 2), Point(0.5, 0.5)], index=[1, 2])
         expected = Series([np.nan, 2.0, np.nan])
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            assert_series_equal(self.g5.project(s), expected)
+        assert_series_equal(self.g5.project(s), expected)
 
         expected = Series([2.0, 0.5], index=self.g5.index)
         assert_series_equal(self.g5.project(s, align=False), expected)
diff --git a/geopandas/tests/test_geoseries.py b/geopandas/tests/test_geoseries.py
index b9f81338..99b8336b 100644
--- a/geopandas/tests/test_geoseries.py
+++ b/geopandas/tests/test_geoseries.py
@@ -110,11 +110,11 @@ class TestSeries:
         # Test that warning is issued when operating on non-aligned series
 
         # _series_op
-        with pytest.warns(UserWarning, match="The indices .+ different"):
+        with pytest.warns(UserWarning, match="The indices .+ not equal"):
             self.a1.contains(self.a2)
 
         # _geo_op
-        with pytest.warns(UserWarning, match="The indices .+ different"):
+        with pytest.warns(UserWarning, match="The indices .+ not equal"):
             self.a1.union(self.a2)
 
     def test_no_warning_if_aligned(self):
@@ -136,8 +136,7 @@ class TestSeries:
         assert_array_equal(self.g1.geom_equals(self.sq), [False, True])
 
     def test_geom_equals_align(self):
-        with pytest.warns(UserWarning, match="The indices .+ different"):
-            a = self.a1.geom_equals(self.a2, align=True)
+        a = self.a1.geom_equals(self.a2, align=True)
         exp = pd.Series([False, True, False], index=["A", "B", "C"])
         assert_series_equal(a, exp)
 
@@ -153,7 +152,7 @@ class TestSeries:
         with warnings.catch_warnings():
             warnings.filterwarnings(
                 "ignore",
-                "The indices of the two GeoSeries are different",
+                "The indices of the left and right GeoSeries' are not equal",
                 UserWarning,
             )
             assert_array_equal(
@@ -170,7 +169,9 @@ class TestSeries:
         assert_array_equal(self.g1.geom_equals_exact(self.sq, 0.001), [False, True])
         with warnings.catch_warnings():
             warnings.filterwarnings(
-                "ignore", "The indices of the two GeoSeries are different", UserWarning
+                "ignore",
+                "The indices of the left and right GeoSeries' are not equal",
+                UserWarning,
             )
             assert_array_equal(
                 self.a1.geom_equals_exact(self.a2, 0.001, align=True),