Compute cost distance measured from each start point to all end points.
        The function returns the distances between the start point and the end
        points as a Pandas dataframe. Additionally, for each start point it computes
        the level of isolation, i.e. its average travel distance to all other locations

def distance(self, sources, destinations, x='x', y='y', isolation=True,
                 export_raster=False, export_shape=False, routes=False, path='./'):
        """
        Compute cost distance measured from each start point to all end points.
        The function returns the distances between the start point and the end
        points as a Pandas dataframe. Additionally, for each start point it computes
        the level of isolation, i.e. its average travel distance to all other locations
        """
        start_points = sources.copy()
        end_points = destinations.copy()
        if (not isinstance(start_points, pd.core.frame.DataFrame) and
            not isinstance(start_points, gp.geodataframe.GeoDataFrame)):
            raise TypeError('Sources has to be a (Geo)Pandas Data Frame Object.')
        if (not isinstance(end_points, pd.core.frame.DataFrame) and
            not isinstance(end_points, gp.geodataframe.GeoDataFrame)):
            raise TypeError('Destinations has to be a (Geo)Pandas Data Frame Object.')
        if not self.mcp_cost:
            self.mcp()
        count = 0
        start_points['row'], start_points['col'] = self.map_pixel_location(start_points[x],
                                                                           start_points[y])
        end_points['row'], end_points['col'] = self.map_pixel_location(end_points[x], end_points[y])
        start_points['ID'] = start_points.index.values
        end_points['ID'] = end_points.index.values+start_points['ID'].max()+1

        for i in start_points.iterrows():
            cumulative_costs, traceback = self.mcp_cost.find_costs([[i[1].row, i[1].col]])
            dist = cumulative_costs[end_points.row.values, end_points.col.values].transpose()/(7*24)
            df2 = pd.DataFrame(np.array([(i[1]['ID']*np.ones_like(dist)).flatten(),
                                    end_points['ID'], dist.flatten()]).transpose(),
                                    columns=['ID1', 'ID2', 'dist'])
            # Keep only locations that are accessible
            df2 = df2.loc[df2['dist'] < np.inf]
            if isolation:
                grisolation = np.ma.masked_array(cumulative_costs,
                                    mask=np.logical_or(self.raster.mask, cumulative_costs == np.inf)
                                    , fill_value=np.nan).mean()/(7*24)
                start_points.loc[i[0], 'Iso'] = grisolation
            if export_raster:
                cumulative_costs = GeoRaster(np.ma.masked_array(cumulative_costs,
                                                                mask=np.logical_or(self.raster.mask,
                                                                                   cumulative_costs == np.inf),
                                                                fill_value=np.nan), self.geot, self.nodata_value,
                                             projection=self.projection, datatype=self.datatype)
                cumulative_costs.raster.data[cumulative_costs.raster.mask] = cumulative_costs.nodata_value
                cumulative_costs.to_tiff(path+str(i[1]['ID']))
            if df2.size > 0:
                if export_shape:
                    routes = True
                if routes:
                    df2['geometry'] = df2['ID2'].apply(lambda x:
                                                self.mcp_cost.traceback(end_points.loc[end_points['ID'] == x][['row', 'col']].values[0]))
                    df2['geometry'] = df2.geometry.apply(lambda x: [map_pixel_inv(y[0], y[1], self.geot[1],
                                                    self.geot[-1], self.geot[0], self.geot[-3]) for y in x])
                    df2['geometry'] = df2.geometry.apply(lambda x: LineString(x) if int(len(x) > 1)
                                                         else LineString([x[0], x[0]]))
                    df2 = gp.GeoDataFrame(df2, crs=cea)
                if isolation:
                    df2['Iso'] = grisolation
                if count == 0:
                    self.grdist = df2.copy()
                else:
                    self.grdist = self.grdist.append(df2)
                count += 1
        if routes:
            self.grdist = gp.GeoDataFrame(self.grdist, crs=cea)
        if export_shape:
            start_pointscols = sources.columns.values
            end_pointscols = destinations.columns.values
            if 'geometry' in end_pointscols:
                self.grdist = pd.merge(self.grdist, end_points[['ID'] + end_pointscols.tolist()].drop('geometry', axis=1), left_on='ID2', right_on='ID', how='left')
            else:
                self.grdist = pd.merge(self.grdist, end_points[['ID']+end_pointscols.tolist()], left_on='ID2', right_on='ID', how='left')
            if 'geometry' in self.start_pointscols:
                self.grdist = pd.merge(self.grdist, start_points[['ID']+start_pointscols.tolist()].drop('geometry', axis=1), left_on='ID1', right_on='ID', how='left',
                             suffixes=['_2', '_1'])
            else:
                self.grdist = pd.merge(self.grdist, start_points[['ID']+start_pointscols.tolist()], left_on='ID1', right_on='ID', how='left',
                             suffixes=['_2', '_1'])
            self.grdist = gp.GeoDataFrame(self.grdist, crs=cea)
            self.grdist.to_file(path+'routes.shp')