diff --git a/docs/read_geotif.md b/docs/read_geotif.md
index 5f9378f..9a05148 100644
--- a/docs/read_geotif.md
+++ b/docs/read_geotif.md
@@ -2,19 +2,20 @@
 
 Read in data (from tif format, most likely georeferenced image data). 
 
-**plantcv.geospatial.read_geotif**(*filename, bands="B,G,R"*)
+**plantcv.geospatial.read_geotif**(*filename, bands="R,G,B", cropto=None*)
 
 **returns** [PlantCV Spectral_data](https://plantcv.readthedocs.io/en/latest/Spectral_data/) object instance.
 
 - **Parameters:**
     - filename - Path of the TIF image file.
-    - bands - Comma separated string representing the order of image bands (e.g., `bands="B,G,R"`), or a list of wavelengths (e.g., `bands=[650, 560, 480]`)
+    - bands - Comma separated string representing the order of image bands (e.g., `bands="R,G,B"`), or a list of wavelengths (e.g., `bands=[650, 560, 480]`)
         - Supported symbols and their default wavelengths: 
             - R (red) = 650nm
             - G (green) = 560nm
             - B (blue) = 480nm
             - RE (rededge) = 717nm
             - N or NIR (near infrared) = 842nm
+    - cropto - A geoJSON-type shapefile to crop the input image as it is read in. Default is None. 
 
 - **Example use:**
     - below
@@ -25,7 +26,8 @@ import plantcv.geospatial as geo
 
 # Read geotif in
 ortho1 = geo.read_geotif(filename="./data/example_img.tif", bands="b,g,r,RE,NIR")
-ortho2 = geo.read_geotif(filename="./data/example_rgb_img.tif", bands="B,G,R")
+ortho2 = geo.read_geotif(filename="./data/example_rgb_img.tif", bands="R,G,B",
+                         cropto="./shapefiles/experimental_bounds.geojson)
 
 ```
 
diff --git a/plantcv/geospatial/read_geotif.py b/plantcv/geospatial/read_geotif.py
index 1e63de4..451aef8 100644
--- a/plantcv/geospatial/read_geotif.py
+++ b/plantcv/geospatial/read_geotif.py
@@ -4,14 +4,16 @@
 import cv2
 import rasterio
 import numpy as np
-from plantcv.plantcv import params
-from plantcv.plantcv import fatal_error
+import fiona
+from rasterio.mask import mask
+from plantcv.plantcv import warn, params, fatal_error
 from plantcv.plantcv._debug import _debug
 from plantcv.plantcv.classes import Spectral_data
+from shapely.geometry import shape, MultiPoint, mapping
 
 
 def _find_closest_unsorted(array, target):
-    """Find closest index of array item with smallest distance from the target.
+    """Find closest index of array item with smallest distance from
 
     Parameters
     ----------
@@ -60,7 +62,7 @@ def _parse_bands(bands):
     return band_list
 
 
-def read_geotif(filename, bands="B,G,R"):
+def read_geotif(filename, bands="R,G,B", cropto=None):
     """Read Georeferenced TIF image from file.
 
     Parameters
@@ -68,30 +70,55 @@ def read_geotif(filename, bands="B,G,R"):
     filename : str
         Path of the TIF image file.
     bands : str, list, optional
-        Comma separated string listing the order of bands or a list of wavelengths, by default "B,G,R"
+        Comma separated string listing the order of bands or a list of wavelengths, by default "R,G,B"
 
     Returns
     -------
     plantcv.plantcv.classes.Spectral_data
-        Orthomosaic image data
+        Orthomosaic image data in a Spectral_data class instance
     """
-    img = rasterio.open(filename)
-    img_data = img.read()
+    if cropto:
+        with fiona.open(cropto, 'r') as shapefile:
+            # polygon-type shapefile
+            if len(shapefile) == 1:
+                shapes = [feature['geometry'] for feature in shapefile]
+            # points-type shapefile
+            else:
+                points = [shape(feature["geometry"]) for feature in shapefile]
+                multi_point = MultiPoint(points)
+                convex_hull = multi_point.convex_hull
+                shapes = [mapping(convex_hull)]
+        # rasterio does the cropping within open
+        with rasterio.open(filename, 'r') as src:
+            img_data, geo_transform = mask(src, shapes, crop=True)
+            d_type = src.dtypes[0]
+            geo_crs = src.crs.wkt
+
+    else:
+        img = rasterio.open(filename)
+        img_data = img.read()
+        d_type = img.dtypes[0]
+        geo_transform = img.transform
+        geo_crs = img.crs.wkt
+
     img_data = img_data.transpose(1, 2, 0)  # reshape such that z-dimension is last
-    height = img.height
-    width = img.width
+    height, width, depth = img_data.shape
     wavelengths = {}
 
     # Parse bands if input is a string
     if isinstance(bands, str):
         bands = _parse_bands(bands)
-
     # Create a dictionary of wavelengths and their indices
     for i, wl in enumerate(bands):
         wavelengths[wl] = i
-
+    # Check if user input matches image dimension in z direction
+    if depth != len(bands):
+        warn(f"{depth} bands found in the image data but {filename} was provided with {bands}")
     # Mask negative background values
     img_data[img_data < 0.] = 0
+    if np.sum(img_data) == 0:
+        fatal_error(f"your image is empty, are the crop-to bounds outside of the {filename} image area?")
+
     # Make a list of wavelength keys
     wl_keys = wavelengths.keys()
     # Find which bands to use for red, green, and blue bands of the pseudo_rgb image
@@ -113,12 +140,15 @@ def read_geotif(filename, bands="B,G,R"):
                                    max_wavelength=None,
                                    min_wavelength=None,
                                    max_value=np.max(img_data), min_value=np.min(img_data),
-                                   d_type=img.dtypes[0],
+                                   d_type=d_type,
                                    wavelength_dict=wavelengths, samples=int(width),
                                    lines=int(height), interleave=None,
                                    wavelength_units="nm", array_type="datacube",
-                                   pseudo_rgb=pseudo_rgb, filename=filename, default_bands=[480, 540, 630])
-
-    _debug(visual=pseudo_rgb, filename=os.path.join(params.debug_outdir, f"{params.device}_pseudo_rgb.png"))
+                                   pseudo_rgb=pseudo_rgb, filename=filename,
+                                   default_bands=[480, 540, 630],
+                                   geo_transform=geo_transform,
+                                   geo_crs=geo_crs)
 
+    _debug(visual=pseudo_rgb, filename=os.path.join(params.debug_outdir,
+                                                    f"{params.device}_pseudo_rgb.png"))
     return spectral_array
diff --git a/pyproject.toml b/pyproject.toml
index b95d666..7617ecb 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -13,7 +13,7 @@ dependencies = [
     "matplotlib",
     "rasterio",
     "fiona",
-
+    "shapely",
 ]
 requires-python = ">=3.6"
 authors = [
diff --git a/tests/conftest.py b/tests/conftest.py
index a1ca38c..c3a7bdf 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -14,10 +14,16 @@ def __init__(self):
         self.snapshot_dir = os.path.join(self.datadir, "snapshot_dir")
         # multispectral image
         self.cropped_tif = os.path.join(self.datadir, "615.tif")
+        # empty image
+        self.empty_tif = os.path.join(self.datadir, "cropped_empty.tif")
         # rgb image
         self.rgb_tif = os.path.join(self.datadir, "rgb.tif")
         # points shapefilex
         self.pts_geojson = os.path.join(self.datadir, "kura_whole_field.shp")
+        # polygon shapefile
+        self.square_crop = os.path.join(self.datadir, "square_crop.geojson")
+        # points shapefile
+        self.point_crop = os.path.join(self.datadir, "point_crop.geojson")
 
 @pytest.fixture(scope="session")
 def test_data():
diff --git a/tests/test_geospatial_read_geotif.py b/tests/test_geospatial_read_geotif.py
index ad1d698..7f268b8 100644
--- a/tests/test_geospatial_read_geotif.py
+++ b/tests/test_geospatial_read_geotif.py
@@ -15,7 +15,7 @@ def test_geospatial_read_geotif_rgb(test_data):
     """Test for plantcv-geospatial."""
     # read in small 5-band tif image
     img = read_geotif(filename=test_data.rgb_tif, bands="R,G,B")
-    assert img.array_data.shape[2] == 4
+    assert img.pseudo_rgb.shape == (284, 261, 3)
 
 
 def test_geospatial_read_geotif_bad_input(test_data):
@@ -23,3 +23,25 @@ def test_geospatial_read_geotif_bad_input(test_data):
     # read in small 5-band tif image
     with pytest.raises(RuntimeError):
         _ = read_geotif(filename=test_data.cropped_tif, bands="p")
+
+
+def test_geospatial_read_geotif_bad_crop(test_data):
+    """Test for plantcv-geospatial."""
+    # read in small 5-band tif image
+    with pytest.raises(RuntimeError):
+        _ = read_geotif(filename=test_data.empty_tif, bands="B,G,R,RE,N")
+
+
+def test_geospatial_read_geotif_polygon_crop(test_data):
+    """Test for plantcv-geospatial."""
+    # read in rgb image with a polygon-type shapefile
+    img = read_geotif(filename=test_data.rgb_tif, bands=[650, 560, 480],
+                      cropto=test_data.square_crop)
+    assert img.pseudo_rgb.shape == (80, 83, 3)
+
+
+def test_geospatial_read_geotif_point_crop(test_data):
+    """Test for plantcv-geospatial."""
+    # read in rgb image with a polygon-type shapefile
+    img = read_geotif(filename=test_data.rgb_tif, bands="R,G,B", cropto=test_data.point_crop)
+    assert img.pseudo_rgb.shape == (41, 46, 3)
diff --git a/tests/testdata/cropped_empty.tif b/tests/testdata/cropped_empty.tif
new file mode 100644
index 0000000..85e00ce
Binary files /dev/null and b/tests/testdata/cropped_empty.tif differ
diff --git a/tests/testdata/point_crop.geojson b/tests/testdata/point_crop.geojson
new file mode 100644
index 0000000..53f469c
--- /dev/null
+++ b/tests/testdata/point_crop.geojson
@@ -0,0 +1,11 @@
+{
+"type": "FeatureCollection",
+"name": "point_crop",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32615" } },
+"features": [
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 720198.797875224030577, 4302928.175875684246421 ] } },
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 720199.770441700122319, 4302928.186563228257 ] } },
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 720199.781129243900068, 4302927.310184645466506 ] } },
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 720198.851312942570075, 4302927.352934820577502 ] } }
+]
+}
diff --git a/tests/testdata/square_crop.geojson b/tests/testdata/square_crop.geojson
new file mode 100644
index 0000000..82b8505
--- /dev/null
+++ b/tests/testdata/square_crop.geojson
@@ -0,0 +1,8 @@
+{
+"type": "FeatureCollection",
+"name": "square_crop",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32615" } },
+"features": [
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 720200.657507826690562, 4302928.218625859357417 ], [ 720199.535315738874488, 4302929.041566723957658 ], [ 720198.855781324207783, 4302928.114928886294365 ], [ 720199.977973412023857, 4302927.291988021694124 ], [ 720200.657507826690562, 4302928.218625859357417 ] ] ] } }
+]
+}