diff --git a/xray/backends/api.py b/xray/backends/api.py
index a6030952d5b..c2201c6ac78 100644
--- a/xray/backends/api.py
+++ b/xray/backends/api.py
@@ -1,6 +1,7 @@
 import sys
 import gzip
 import itertools
+import re
 from glob import glob
 from io import BytesIO
 
@@ -14,17 +15,29 @@
 from ..core.pycompat import basestring, OrderedDict, range
 
 
-def _get_default_netcdf_engine(engine):
-    try:
-        import netCDF4
-        engine = 'netcdf4'
-    except ImportError: # pragma: no cover
+def _get_default_engine(path, allow_remote=False):
+    if allow_remote and re.search('^https?\://', path):  # pragma: no cover
         try:
-            import scipy.io.netcdf
-            engine = 'scipy'
+            import netCDF4
+            engine = 'netcdf4'
         except ImportError:
-            raise ValueError('cannot read or write netCDF files without '
-                             'netCDF4-python or scipy installed')
+            try:
+                import pydap
+                engine = 'pydap'
+            except ImportError:
+                raise ValueError('netCDF4 or pydap is required for accessing '
+                                 'remote datasets via OPeNDAP')
+    else:
+        try:
+            import netCDF4
+            engine = 'netcdf4'
+        except ImportError: # pragma: no cover
+            try:
+                import scipy.io.netcdf
+                engine = 'scipy'
+            except ImportError:
+                raise ValueError('cannot read or write netCDF files without '
+                                 'netCDF4-python or scipy installed')
     return engine
 
 
@@ -67,7 +80,7 @@ def open_dataset(filename_or_obj, group=None, decode_cf=True,
     engine : {'netcdf4', 'scipy', 'pydap', 'h5netcdf'}, optional
         Engine to use when reading netCDF files. If not provided, the default
         engine is chosen based on available dependencies, with a preference for
-        'netcdf4' if reading a file on disk.
+        'netcdf4'.
     blockdims, blockshape : dict, optional
         If blockdims or blockshape is provided, it used to load the new dataset
         into dask arrays. This is an experimental feature; see the
@@ -116,7 +129,8 @@ def maybe_decode_store(store):
             # TODO: automatically fall back to using pydap if given a URL and
             # netCDF4 is not available
             if engine is None:
-                engine = _get_default_netcdf_engine(engine)
+                engine = _get_default_engine(filename_or_obj,
+                                             allow_remote=True)
             if engine == 'netcdf4':
                 store = backends.NetCDF4DataStore(filename_or_obj, group=group)
             elif engine == 'scipy':
@@ -198,7 +212,7 @@ def to_netcdf(dataset, path=None, mode='w', format=None, group=None,
                              'to_netcdf: %r. Only the default engine '
                              "or engine='scipy' is supported" % engine)
     elif engine is None:
-        engine = _get_default_netcdf_engine(engine)
+        engine = _get_default_engine(path)
 
     write_funcs = {'netcdf4': _to_netcdf4,
                    'scipy': _to_scipy_netcdf,
diff --git a/xray/backends/pydap_.py b/xray/backends/pydap_.py
index 49dda7ff065..d741e0722ea 100644
--- a/xray/backends/pydap_.py
+++ b/xray/backends/pydap_.py
@@ -26,14 +26,17 @@ def __getitem__(self, key):
         if not isinstance(key, tuple):
             key = (key,)
         for k in key:
-            if not (isinstance(k, int)
-                    or isinstance(k, slice)
-                    or k is Ellipsis):
+            if not (isinstance(k, (int, np.integer, slice)) or k is Ellipsis):
                 raise IndexError('pydap only supports indexing with int, '
                                  'slice and Ellipsis objects')
         # pull the data from the array attribute if possible, to avoid
         # downloading coordinate data twice
-        return getattr(self.array, 'array', self.array)[key]
+        array = getattr(self.array, 'array', self.array)
+        result = array[key]
+        # pydap doesn't squeeze axes automatically like numpy
+        axis = tuple(k for k in key if isinstance(k, (int, np.integer)))
+        result = np.squeeze(result, axis=axis)
+        return result
 
 
 class PydapDataStore(AbstractDataStore):
diff --git a/xray/test/test_backends.py b/xray/test/test_backends.py
index 365d97dbb38..c9bf064819d 100644
--- a/xray/test/test_backends.py
+++ b/xray/test/test_backends.py
@@ -714,15 +714,29 @@ def test_dask_roundtrip(self):
                         self.assertDatasetIdentical(data, on_disk)
 
 
-@requires_netCDF4
+@requires_scipy_or_netCDF4
 @requires_pydap
 class PydapTest(TestCase):
     def test_cmp_local_file(self):
         url = 'http://test.opendap.org/opendap/hyrax/data/nc/bears.nc'
-        actual = open_dataset(url, engine='pydap')
-        with open_example_dataset('bears.nc') as expected:
-            # don't check attributes since pydap doesn't serialize them correctly
-            # also skip the "bears" variable since the test DAP server incorrectly
-            # concatenates it.
-            self.assertDatasetEqual(actual.drop('bears'),
-                                    expected.drop('bears'))
+
+        @contextlib.contextmanager
+        def create_datasets():
+            actual = open_dataset(url, engine='pydap')
+            with open_example_dataset('bears.nc') as expected:
+                # don't check attributes since pydap doesn't serialize them
+                # correctly also skip the "bears" variable since the test DAP
+                # server incorrectly concatenates it.
+                actual = actual.drop('bears')
+                expected = expected.drop('bears')
+                yield actual, expected
+
+        with create_datasets() as (actual, expected):
+            self.assertDatasetEqual(actual, expected)
+
+        with create_datasets() as (actual, expected):
+            self.assertDatasetEqual(actual.isel(i=0), expected.isel(i=0))
+
+        with create_datasets() as (actual, expected):
+            self.assertDatasetEqual(actual.isel(j=slice(1, 2)),
+                                    expected.isel(j=slice(1, 2)))