diff --git a/ci/requirements/py36-hypothesis.yml b/ci/requirements/py36-hypothesis.yml
index f1aad7b39e0..57f4a236845 100644
--- a/ci/requirements/py36-hypothesis.yml
+++ b/ci/requirements/py36-hypothesis.yml
@@ -22,6 +22,7 @@ dependencies:
   - seaborn
   - toolz
   - rasterio
+  - boto3
   - bottleneck
   - zarr
   - pydap
diff --git a/ci/requirements/py36.yml b/ci/requirements/py36.yml
index a8481623f80..60259d59965 100644
--- a/ci/requirements/py36.yml
+++ b/ci/requirements/py36.yml
@@ -22,6 +22,7 @@ dependencies:
   - seaborn
   - toolz
   - rasterio
+  - boto3
   - bottleneck
   - zarr
   - pseudonetcdf>=3.0.1
diff --git a/ci/requirements/py37-windows.yml b/ci/requirements/py37-windows.yml
index 748546c75b3..51936279134 100644
--- a/ci/requirements/py37-windows.yml
+++ b/ci/requirements/py37-windows.yml
@@ -22,4 +22,5 @@ dependencies:
   - seaborn
   - toolz
   - rasterio
+  - boto3
   - zarr
diff --git a/ci/requirements/py37.yml b/ci/requirements/py37.yml
index de7fa3abdc6..f1b6d46fd95 100644
--- a/ci/requirements/py37.yml
+++ b/ci/requirements/py37.yml
@@ -23,6 +23,7 @@ dependencies:
   - seaborn
   - toolz
   - rasterio
+  - boto3
   - bottleneck
   - zarr
   - pseudonetcdf>=3.0.1
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index 2235be4fbb7..332d1b3b15a 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -3743,7 +3743,6 @@ def test_rasterio_vrt(self):
                 with rasterio.vrt.WarpedVRT(src, crs='epsg:4326') as vrt:
                     expected_shape = (vrt.width, vrt.height)
                     expected_crs = vrt.crs
-                    print(expected_crs)
                     expected_res = vrt.res
                     # Value of single pixel in center of image
                     lon, lat = vrt.xy(vrt.width // 2, vrt.height // 2)
@@ -3751,7 +3750,6 @@ def test_rasterio_vrt(self):
                     with xr.open_rasterio(vrt) as da:
                         actual_shape = (da.sizes['x'], da.sizes['y'])
                         actual_crs = da.crs
-                        print(actual_crs)
                         actual_res = da.res
                         actual_val = da.sel(dict(x=lon, y=lat),
                                             method='nearest').data
@@ -3790,35 +3788,29 @@ def test_rasterio_vrt_with_transform_and_size(self):
 
     @network
     def test_rasterio_vrt_network(self):
+        # Make sure loading w/ rasterio give same results as xarray
         import rasterio
-
-        url = 'https://storage.googleapis.com/\
-        gcp-public-data-landsat/LC08/01/047/027/\
-        LC08_L1TP_047027_20130421_20170310_01_T1/\
-        LC08_L1TP_047027_20130421_20170310_01_T1_B4.TIF'
-        env = rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN='EMPTY_DIR',
-                           CPL_VSIL_CURL_USE_HEAD=False,
-                           CPL_VSIL_CURL_ALLOWED_EXTENSIONS='TIF')
-        with env:
-            with rasterio.open(url) as src:
+        # use same url that rasterio package uses in tests
+        prefix = "https://landsat-pds.s3.amazonaws.com/L8/139/045/"
+        image = "LC81390452014295LGN00/LC81390452014295LGN00_B1.TIF"
+        httpstif = prefix + image
+        with rasterio.Env(aws_unsigned=True):
+            with rasterio.open(httpstif) as src:
                 with rasterio.vrt.WarpedVRT(src, crs='epsg:4326') as vrt:
-                    expected_shape = (vrt.width, vrt.height)
-                    expected_crs = vrt.crs
+                    expected_shape = vrt.width, vrt.height
                     expected_res = vrt.res
                     # Value of single pixel in center of image
                     lon, lat = vrt.xy(vrt.width // 2, vrt.height // 2)
                     expected_val = next(vrt.sample([(lon, lat)]))
                     with xr.open_rasterio(vrt) as da:
-                        actual_shape = (da.sizes['x'], da.sizes['y'])
-                        actual_crs = da.crs
+                        actual_shape = da.sizes['x'], da.sizes['y']
                         actual_res = da.res
                         actual_val = da.sel(dict(x=lon, y=lat),
                                             method='nearest').data
 
-                        assert_equal(actual_shape, expected_shape)
-                        assert_equal(actual_crs, expected_crs)
-                        assert_equal(actual_res, expected_res)
-                        assert_equal(expected_val, actual_val)
+                        assert actual_shape == expected_shape
+                        assert actual_res == expected_res
+                        assert expected_val == actual_val
 
 
 class TestEncodingInvalid: