diff --git a/.gitignore b/.gitignore
index 102eb729..c99095bc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -161,6 +161,7 @@ cython_debug/
 #.idea/
 .vscode/*
 tests/.vscode/*
+test/
 dnora/.vscode/*
 tests/dnora_wind_temp/*
 tests/dnora_current_temp/*
@@ -177,3 +178,4 @@ grd/
 dnplot/
 spc/
 wlv/
+dnora_bnd_temp/
\ No newline at end of file
diff --git a/dnora/file_module.py b/dnora/file_module.py
index 23e42130..634194cd 100644
--- a/dnora/file_module.py
+++ b/dnora/file_module.py
@@ -352,7 +352,7 @@ def replace_objects(
     """
     for obj_type, obj_name in dict_of_object_names.items():
         if obj_name is not None:
-            filename = re.sub(f"#{obj_type.name}", obj_name, filename, 1)
+            filename = re.sub(f"#{obj_type.name}", obj_name, filename, count=1)
 
     return filename
 
diff --git a/dnora/modelrun/import_functions.py b/dnora/modelrun/import_functions.py
index 062235c4..913e45a1 100644
--- a/dnora/modelrun/import_functions.py
+++ b/dnora/modelrun/import_functions.py
@@ -12,7 +12,11 @@
 from dnora.pick import PointPicker, Area, NearestGridPoint
 from geo_skeletons import PointSkeleton
 from dnora import msg
+from dnora.utils.distance import clustered_around_lon180, wrapped_lon_edges
+from dnora.utils.grid import cluster_points
 
+from dnora.grid import TriGrid
+from dnora.grid.mask import All
 
 def import_data(
     grid: Grid,
@@ -28,6 +32,7 @@ def import_data(
     filename: str,
     point_mask=None,
     point_picker=None,
+    max_calls: int = None,
     **kwargs,
 ) -> DnoraObject:
     """Imports data using DataReader and creates and returns a DNORA object"""
@@ -35,51 +40,130 @@ def import_data(
     msg.print_line(marker="+")
     msg.plain(f"Starting import of data: {obj_type.name}")
     msg.print_line(marker="+")
-    msg.plain(
-        f"Area: {grid.core.x_str}: {grid.edges('lon',native=True)}, {grid.core.y_str}: {grid.edges('lat',native=True)}"
+
+    if grid.core.x_str == 'lon' and clustered_around_lon180(grid.lon()):
+        lon0 = float(np.min(grid.lon()[grid.lon() > 0]))
+        lon1 = float(np.max(grid.lon()[grid.lon() < 0]))
+        msg.plain(
+        f"Area: {grid.core.x_str}: {(lon0, lon1)}, "
+        f"{grid.core.y_str}: {grid.edges('lat',native=True)}"
     )
+    else:
+        msg.plain(
+            f"Area: {grid.core.x_str}: {grid.edges('lon',native=True)}, "
+            f"{grid.core.y_str}: {grid.edges('lat',native=True)}"
+        )
+    
+    if max_calls is not None:
+        msg.plain(f"Max calls set to {max_calls}. Clustering will be used if needed.")
+        clustered_points = cluster_points(lon = grid.lon(), lat = grid.lat(), N_cluster = max_calls)
+        msg.plain(f"Number of clusters created: {len(clustered_points)}")
+        for i, cluster in enumerate(clustered_points):
+            if clustered_around_lon180(cluster[:,0]):
+                lon0 = float(np.min(cluster[cluster[:,0] > 0][:,0]))
+                lon1 = float(np.max(cluster[cluster[:,0] < 0][:,0]))
+                msg.plain(f"Area cluster {i+1} (wrapped): lon: {lon0:.2f} to {lon1:.2f}, lat: {cluster[:,1].min():.2f} to {cluster[:,1].max():.2f}")
+            else:   
+                msg.plain(f"Area cluster {i+1}: lon: {cluster[:,0].min():.2f} to {cluster[:,0].max():.2f}, lat: {cluster[:,1].min():.2f} to {cluster[:,1].max():.2f}")
+
     msg.plain(f"{start_time} - {end_time}")
 
     if dry_run:
         msg.info("Dry run! No data will be imported.")
         return
+    if max_calls is not None:
+        if dnora_objects.get(obj_type).is_gridded():
+            raise ValueError("Cannot use 'max_calls' when importing gridded objects!")
 
-    if not dnora_objects.get(obj_type).is_gridded():
         msg.header(point_picker, "Choosing points to import...")
-        inds = pick_points(
-            grid,
+        msg.plain("Point picking with clustering...")
+
+        obj_list = []
+        clustered_points = [clustered_points[-1]]
+        for i, cluster in enumerate(clustered_points):
+            msg.plain(f"\nPicking points for cluster {i+1}...")
+            lon,lat = cluster[:,0], cluster[:,1]
+            
+            points = TriGrid(lon=lon[:], lat=lat[:], name=f"{grid.name}_cluster_{i}")
+            points.set_boundary_points(All())
+            inds = pick_points(
+                points,
+                reader,
+                start_time,
+                point_picker,
+                expansion_factor,
+                points.boundary_mask(),
+                source,
+                folder,
+                filename,
+                **kwargs,
+            )
+            if len(inds) < 1:
+                msg.warning("PointPicker didn't find any points. Aborting import of data.")
+                return
+
+            msg.plain(f"\nImporting data for cluster {i+1}...")
+            obj = read_data_and_create_object(
+                obj_type,
+                reader,
+                expansion_factor,
+                points,
+                start_time,
+                end_time,
+                name,
+                source,
+                folder,
+                filename,
+                inds,
+                **kwargs,
+            )
+            obj_list.append(obj)
+
+        msg.plain("Concatenating data from all clusters...")
+        obj = None
+        for ob in obj_list:
+            if obj is None:
+                obj = ob
+            else:
+                obj = obj.absorb(ob, dim='inds')
+        msg.plain("Concatenation done.")
+
+
+    else:
+        if not dnora_objects.get(obj_type).is_gridded():
+            msg.header(point_picker, "Choosing points to import...")
+            inds = pick_points(
+                grid,
+                reader,
+                start_time,
+                point_picker,
+                expansion_factor,
+                point_mask,
+                source,
+                folder,
+                filename,
+                **kwargs,
+            )
+            if len(inds) < 1:
+                msg.warning("PointPicker didn't find any points. Aborting import of data.")
+                return
+        else:
+            inds = np.array([])
+        msg.header(reader, f"Importing {obj_type.name}...")
+        obj = read_data_and_create_object(
+            obj_type,
             reader,
-            start_time,
-            point_picker,
             expansion_factor,
-            point_mask,
+            grid,
+            start_time,
+            end_time,
+            name,
             source,
             folder,
             filename,
+            inds,
             **kwargs,
         )
-        if len(inds) < 1:
-            msg.warning("PointPicker didn't find any points. Aborting import of data.")
-            return
-    else:
-        inds = np.array([])
-
-    msg.header(reader, f"Importing {obj_type.name}...")
-
-    obj = read_data_and_create_object(
-        obj_type,
-        reader,
-        expansion_factor,
-        grid,
-        start_time,
-        end_time,
-        name,
-        source,
-        folder,
-        filename,
-        inds,
-        **kwargs,
-    )
     return obj
 
 
@@ -193,31 +277,29 @@ def pick_points(
         y=available_points.get("y"),
     )
 
+    if not np.all(np.logical_not(point_mask)):
+        interest_points = PointSkeleton.from_skeleton(grid, mask=point_mask)
+        slat = interest_points.edges("lat")
+        slon = wrapped_lon_edges(interest_points)
+    else:
+        interest_points = None
+        slat = grid.edges("lat")
+        slon = wrapped_lon_edges(grid)
     ## Only take points that are reasonable close to the wanted grid
     ## This speeds up the searcg considerably, especially if we have points over 84 lat
     ## since then the fast cartesian searhc is not possible
     ## Set a limit for angles close to 90 and -90 lat and -180 and 180 longitude
-    slon, slat = grid.edges("lon"), grid.edges("lat")
-    eps = 1e-12
-    search_grid = Grid(lat=(max(slat[0] - 3, -90 + eps), min(slat[1] + 3, 90 - eps)),
-                        lon=(max(slon[0] - 6, -180 + eps), min(slon[1] + 6, 180 - eps)))
     if isinstance(point_picker, NearestGridPoint):
-        search_inds = Area()(search_grid, all_points, expansion_factor=1)
+        eps = 1e-12
+        slon = (slon[0] - 6, slon[1] + 6)
+
+        search_grid = Grid(lat=(max(slat[0] - 3, -90 + eps), min(slat[1] + 3, 90 - eps)),
+                        lon=(max(slon[0] - 6, -180 + eps), min(slon[1] + 6, 180 - eps))) 
+        search_inds = Area()(search_grid, all_points, expansion_factor=1, lon=slon)
     else:
         search_inds = all_points.inds()
 
-    # if np.all(np.logical_not(point_mask)):
-    #     msg.warning(
-    #         "None of the points set to interest points! Aborting import of data."
-    #     )
-    #     return
-    # else:
-    #     interest_points = PointSkeleton.from_skeleton(grid, mask=point_mask)
 
-    if not np.all(np.logical_not(point_mask)):
-        interest_points = PointSkeleton.from_skeleton(grid, mask=point_mask)
-    else:
-        interest_points = None
 
     inds = point_picker(
         grid=grid,
@@ -227,7 +309,7 @@ def pick_points(
         fast=True,
         **kwargs,
     )
-
     if len(inds) < 1:
         return np.array([])
+
     return search_inds[inds]
diff --git a/dnora/modelrun/modelrun.py b/dnora/modelrun/modelrun.py
index e86d0c41..fa6c7133 100644
--- a/dnora/modelrun/modelrun.py
+++ b/dnora/modelrun/modelrun.py
@@ -19,7 +19,7 @@
 from dnora.type_manager.model_formats import ModelFormat
 from dnora.spectral_grid import SpectralGrid
 
-import dnplot
+
 from dnora import msg
 from dnora.cacher.cache_decorator import cached_reader
 
@@ -159,7 +159,13 @@ def __init__(
         self._nest = {}
         self._parent = None
 
-        self.plot = dnplot.Matplotlib(self)
+        try:
+            import dnplot
+            self.plot = dnplot.Matplotlib(self)
+        except ImportError:
+            msg.info(
+                "dnplot is not installed. Please install dnplot to enable plotting functionality.")
+            
         self._dnora_objects: dict[DnoraDataType, DnoraObject] = {
             DnoraDataType.GRID: grid,
         }
@@ -295,6 +301,7 @@ def _import_data(
         point_mask=None,
         point_picker=None,
         post_process: bool = True,
+        max_calls: Optional[int] = None,
         **kwargs,
     ):
         """Performs import and returns DNORA object"""
@@ -370,6 +377,7 @@ def _import_data(
             filename=filename_to_use,
             point_picker=point_picker,
             point_mask=point_mask,
+            max_calls=max_calls,
             **kwargs,
         )
 
@@ -478,6 +486,7 @@ def import_spectra(
         source: Union[str, DataSource] = DataSource.UNDEFINED,
         folder: Optional[str] = None,
         filename: Optional[str] = None,
+        max_calls: Optional[int] = None,
         **kwargs,
     ) -> None:
         self._import_data(
@@ -491,6 +500,7 @@ def import_spectra(
             filename,
             point_mask=self.grid().boundary_mask(),
             point_picker=point_picker,
+            max_calls=max_calls,
             **kwargs,
         )
 
diff --git a/dnora/pick/point_pickers.py b/dnora/pick/point_pickers.py
index f0bb06cd..c5d4eb86 100644
--- a/dnora/pick/point_pickers.py
+++ b/dnora/pick/point_pickers.py
@@ -3,7 +3,7 @@
 from abc import ABC, abstractmethod
 from typing import Union
 from geo_skeletons import PointSkeleton
-
+from dnora.utils.distance import clustered_around_lon180
 from typing import TYPE_CHECKING
 
 if TYPE_CHECKING:
@@ -93,11 +93,12 @@ def __call__(
         grid: Union[Grid, TriGrid],
         all_points: PointSkeleton,
         expansion_factor: float = 1.5,
+        lon: tuple = None, 
         **kwargs,
     ) -> np.ndarray:
 
         msg.process(f"Using expansion_factor = {expansion_factor:.2f}")
-
+        
         # Define area to search in
         if grid.core.is_cartesian():
             number, zone = grid.utm()
@@ -106,16 +107,30 @@ def __call__(
                 grid.edges("x"), grid.edges("y"), expansion_factor
             )
             x_all, y_all = all_points.xy()
+            maskx = np.logical_and(x_all >= x[0], x_all <= x[1]) 
+    
         else:
+            # If we want to pass a search grid that goes beyond -180 180
+            lat = grid.edges("lat")
+            lon = lon or grid.edges("lon")
             x, y = utils.grid.expand_area(
-                grid.edges("lon"), grid.edges("lat"), expansion_factor
+                lon, lat, expansion_factor, cross_180=True
             )
             x_all, y_all = all_points.lonlat()
 
-        maskx = np.logical_and(x_all >= x[0], x_all <= x[1])
+            if x[0]>x[1]:
+                maskx = np.logical_or(x_all <= x[0], x_all >= x[1]) 
+            else:
+                maskx = np.logical_and(x_all >= x[0], x_all <= x[1]) 
+                if x[0] < -180:
+                    maskx = np.logical_or(maskx, x_all > x[0] + 360)
+                if x[1] > 180:
+                    maskx = np.logical_or(maskx, x_all < x[1] - 360)
+
+
         masky = np.logical_and(y_all >= y[0], y_all <= y[1])
         mask = np.logical_and(maskx, masky)
-
+        
         inds = np.where(mask)[0]
 
         msg.info(
diff --git a/dnora/read/spectra/ec.py b/dnora/read/spectra/ec.py
index 0b8b5dc2..80ec3037 100644
--- a/dnora/read/spectra/ec.py
+++ b/dnora/read/spectra/ec.py
@@ -5,6 +5,7 @@
 from typing import Tuple
 import pandas as pd
 from dnora.process.spectra import RemoveEmpty
+from dnora.utils.distance import assert_lon_almost_equal, clustered_around_lon180
 
 # Import abstract classes and needed instances of them
 from dnora.read.abstract_readers import SpectralDataReader
@@ -122,16 +123,29 @@ def get_coordinates(
         self, grid, start_time, source: DataSource, folder: str, **kwargs
     ) -> dict:
         """Reads first time instance of first file to get longitudes and latitudes for the PointPicker"""
+
+        eps = 1e-12
         lon = np.floor(np.array(grid.edges("lon")) / self.dlon) * self.dlon
         lat = np.floor(np.array(grid.edges("lat")) / self.dlat) * self.dlat
 
-        self._given_grid = Grid(
-            lon=(lon[0] - self.dlon, lon[-1] + self.dlon),
-            lat=(lat[0] - self.dlat, lat[-1] + self.dlat),
-        )
-        self._given_grid.set_spacing(dlon=self.dlon, dlat=self.dlat)
+        if checklon := clustered_around_lon180(grid.lon()): #if the longitudes are centered around 180
+            lon = grid.lon()
+            lon = lon % 360 - 180
+            lon = np.floor(np.array([np.min(lon),np.max(lon)]) / self.dlon) * self.dlon
+            self._given_grid = Grid(lon= (lon[0] - self.dlon, lon[1] + self.dlon),
+                                    lat= (min(max(lat[0] - self.dlat, -90 + eps), 90 - eps),
+                                    min(max(lat[1] + self.dlat, -90 + eps), 90 - eps)))
+        else:        
+            self._given_grid = Grid(
+                lon=(min(max(lon[0] - self.dlon, -180 + eps), 180 - eps),
+                    min(max(lon[-1] + self.dlon, -180 + eps), 180 -eps)),
+                lat=(min(max(lat[0] - self.dlat, -90 + eps), 90 - eps),
+                    min(max(lat[-1] + self.dlat, -90 + eps), 90 - eps)) 
+            )
+        self._given_grid.set_spacing(dlon=self.dlon, dlat=self.dlat) 
         lon_all, lat_all = self._given_grid.lonlat()
-
+        if checklon:
+            lon_all = lon_all % 360 - 180
         return {"lat": lat_all, "lon": lon_all}
 
     def __call__(
@@ -156,16 +170,33 @@ def __call__(
             msg.plain("Removing old files from temporary folder...")
             for f in glob.glob(f"{temp_folder}/EC_ERA5.nc"):
                 os.remove(f)
-
-            nc_file = download_era5_from_cds(
-                start_time,
-                end_time,
-                lon=self._given_grid.edges("lon"),
-                lat=self._given_grid.edges("lat"),
-                dlon=self.dlon,
-                dlat=self.dlat,
-                folder=temp_folder,
-            )
+            if checklon := clustered_around_lon180(grid.lon()): #if the longitudes are centered around 180
+                dlon = self.dlon
+                lon0 = np.min(grid.lon()[grid.lon() > 0])
+                lon1 = np.max(grid.lon()[grid.lon() < 0])
+                
+                lon0 = np.floor(lon0 / dlon) * dlon - dlon
+                lon1 = np.floor(lon1 / dlon) * dlon + dlon
+
+                nc_file = download_era5_from_cds(
+                    start_time,
+                    end_time,
+                    lon=(lon0, lon1),
+                    lat=self._given_grid.edges("lat"),
+                    dlon=self.dlon,
+                    dlat=self.dlat,
+                    folder=temp_folder,
+                )
+            else:
+                nc_file = download_era5_from_cds(
+                    start_time,
+                    end_time,
+                    lon=self._given_grid.edges("lon"),
+                    lat=self._given_grid.edges("lat"),
+                    dlon=self.dlon,
+                    dlat=self.dlat,
+                    folder=temp_folder,
+                )
         else:
             nc_file = f"{temp_folder}/EC_ERA5.nc"
         bnd_spec = xr.open_dataset(nc_file)
@@ -199,7 +230,10 @@ def __call__(
 
         # Check that the points decoded from the file is consistent with the created points given to the PointPicker
         glon, glat = self._given_grid.lonlat()
-        np.testing.assert_array_almost_equal(lon, glon)
+        if checklon:
+            glon = glon % 360 - 180
+            lon = (lon + 180) % 360 -180
+        assert_lon_almost_equal(lon, glon)
         np.testing.assert_array_almost_equal(lat, glat)
 
         # Inds given by point picker
diff --git a/dnora/utils/distance.py b/dnora/utils/distance.py
index 3fc51b77..5aafa835 100644
--- a/dnora/utils/distance.py
+++ b/dnora/utils/distance.py
@@ -47,3 +47,38 @@ def domain_size_in_km(
     km_y = distance_2points(lat[0], lon[0], lat[1], lon[0])
 
     return km_x, km_y
+
+def assert_lon_almost_equal(lon1, lon2, decimal=6):
+    """
+    Assert that two longitude arrays are almost equal,
+    accounting for wraparound at ±180°.
+    """
+    lon1 = np.asarray(lon1)
+    lon2 = np.asarray(lon2)
+
+    diff = (lon1 - lon2 + 180) % 360 - 180
+    np.testing.assert_array_almost_equal(diff, 0.0, decimal=decimal)
+
+def clustered_around_lon180(lon: np.ndarray) -> bool:
+    """
+    Check if longitudes are clustered around the ±180° meridian.
+    """
+    lon = np.asarray(lon)
+
+    neg = lon[lon < 0]
+    pos = lon[lon >= 0]
+
+    # Must have values on both sides
+    if neg.size == 0 or pos.size == 0:
+        return False
+
+    return np.all(neg < -90) and np.all(pos > 90)
+
+def wrapped_lon_edges(grid) -> tuple:
+    """Returns normal edges for normal grids and edges wrapped around +-180 if necessary"""
+    if not clustered_around_lon180(grid.lon()):
+        return grid.edges('lon')
+    else:
+        lon0 = float(np.min(grid.lon()[grid.lon() > 0]))
+        lon1 = float(np.max(grid.lon()[grid.lon() < 0]))
+        return (lon0, lon1)
\ No newline at end of file
diff --git a/dnora/utils/grid.py b/dnora/utils/grid.py
index 91133a15..5a9b63d9 100644
--- a/dnora/utils/grid.py
+++ b/dnora/utils/grid.py
@@ -1,6 +1,7 @@
 import numpy as np
-
-
+from sklearn.cluster import KMeans
+import pandas as pd
+from typing import Union
 def data_covers_grid(skeleton, grid):
     """Checks if a given skeleton covers a given grid"""
     for coord in ["lon", "lat"]:
@@ -142,6 +143,7 @@ def expand_area(
     expansion_factor: float,
     dlon: float = 0.0,
     dlat: float = 0.0,
+    cross_180: bool = False,
 ) -> tuple[float, float, float, float]:
     """
     Expands a lon-lat bounding box with an expansion factor.
@@ -153,8 +155,9 @@ def expand_area(
     expansion_factor=1.2 gives (59.9, 61.1)
     expansion_factor=1.2 and dlat = 0.25 gives (59.75, 61.25)
     """
-
-    expand_lon = (lon[1] - lon[0]) * (expansion_factor - 1) * 0.5
+    #delta_lon = ((lon[1] - lon[0] + 180) % 360) - 180
+    delta_lon = lon[1] - lon[0]
+    expand_lon = (delta_lon) * (expansion_factor - 1) * 0.5
     expand_lat = (lat[1] - lat[0]) * (expansion_factor - 1) * 0.5
 
     expand_lon = np.maximum(expand_lon, dlon)
@@ -162,7 +165,9 @@ def expand_area(
 
     new_lon = lon[0] - expand_lon, lon[1] + expand_lon
     new_lat = lat[0] - expand_lat, lat[1] + expand_lat
-
+    if not cross_180:
+        new_lon = max(new_lon[0], -180.0), min(new_lon[1], 180.0)
+    new_lat = max(new_lat[0], -90.0), min(new_lat[1], 90.0)
     return new_lon, new_lat
 
 
@@ -199,3 +204,33 @@ def get_coordinates_from_ds(ds, return_dict: bool = False) -> tuple:
 
 def all_none(val) -> bool:
     return not [a for a in val if a is not None]
+
+
+def cluster_points(lon: np.ndarray,
+                   lat: np.ndarray,
+                   N_cluster) -> list[np.ndarray]:
+    """Clusters boundary points into N clusters using KMeans.
+    Returns a list of arrays, each array containing the lon-lat points of a cluster.
+    """
+    points = np.column_stack((lon, lat))
+    
+    # Convert lon/lat → XYZ unit sphere
+    lon_rad = np.radians(points[:,0])
+    lat_rad = np.radians(points[:,1])
+
+    x = np.cos(lat_rad) * np.cos(lon_rad)
+    y = np.cos(lat_rad) * np.sin(lon_rad)
+    z = np.sin(lat_rad)
+
+    points_xyz = np.column_stack((x, y, z))
+
+    # ---- KMEANS ON 3D XYZ ----
+    kmeans = KMeans(n_clusters=N_cluster, random_state=0).fit(points_xyz)
+    labels = kmeans.labels_
+
+    clustered_points = []
+    for cluster_id in range(N_cluster):
+        cluster_points = points[labels == cluster_id]
+        clustered_points.append(cluster_points)
+    return clustered_points
+
diff --git a/dnora/utils/io.py b/dnora/utils/io.py
index df6bf739..6158200b 100644
--- a/dnora/utils/io.py
+++ b/dnora/utils/io.py
@@ -64,9 +64,9 @@ def get_url(
             else:
                 url_temp = str(url_temp)
 
-            url_temp = re.sub("https:/", "https://", url_temp, 1)
-            url_temp = re.sub("http:/", "http://", url_temp, 1)
-            url_temp = re.sub("ftp:/", "ftp://", url_temp, 1)
+            url_temp = re.sub("https:/", "https://", url_temp, count=1)
+            url_temp = re.sub("http:/", "http://", url_temp, count=1)
+            url_temp = re.sub("ftp:/", "ftp://", url_temp, count=1)
             if time_stamp is not None:
                 for floor_hour in range(1, 24):
                     hfloor = int(np.floor(time_stamp.hour / floor_hour) * floor_hour)
diff --git a/dnora/utils/spec.py b/dnora/utils/spec.py
index ca206e95..1dff118e 100644
--- a/dnora/utils/spec.py
+++ b/dnora/utils/spec.py
@@ -1,7 +1,7 @@
 import numpy as np
 from scipy import interpolate
 import scipy
-
+import xarray as xr
 
 def directional_distribution(freq, fp, dirs, dirp):
     """Calculates directional cos**2s(0.5*theta) distribution of spectrum
diff --git a/dnora/utils/time.py b/dnora/utils/time.py
index 348bb140..90a38c4a 100644
--- a/dnora/utils/time.py
+++ b/dnora/utils/time.py
@@ -141,20 +141,28 @@ def create_time_stamps(
     end_times = start_times + pd.DateOffset(hours=stride - 1)
 
     # First time might not coincide with first step in first file
-    start_times.values[0] = pd.Timestamp(start_time)
-
+    
+    #start_times.values[0] = pd.Timestamp(start_time)
+    start_times = start_times.insert(0, pd.Timestamp(start_time)).delete(1)
     if last_file:
         # In operational systems we might want to read a longer segment from the last file
-        end_times.values[-1] = min(
+        new_end_time = min(
             [
                 pd.Timestamp(last_file) + pd.DateOffset(hours=(hours_per_file - 1)),
                 pd.Timestamp(end_time),
             ]
         )
     else:
-        # Last time might not coincide with last step in last file
-        end_times.values[-1] = pd.Timestamp(end_time)
-
+        new_end_time = pd.Timestamp(end_time)
+        
+        # end_times.values[-1] = min(
+        #     [
+        #         pd.Timestamp(last_file) + pd.DateOffset(hours=(hours_per_file - 1)),
+        #         pd.Timestamp(end_time),
+        #     ]
+        # )
+
+    end_times = end_times[:-1].append(pd.DatetimeIndex([new_end_time]))
     return start_times, end_times, file_times
 
 def get_first_file(
diff --git a/environment.yml b/environment.yml
index 84c04d1d..eeff80ae 100644
--- a/environment.yml
+++ b/environment.yml
@@ -10,7 +10,7 @@ dependencies:
   - xarray>=0.20.1
   - netcdf4
   - meshio
-  - geo-skeletons>=0.21.1
+  - geo-skeletons=0.21.1
   - geo-parameters>=0.11.0
   - pytest
   - pip
@@ -19,3 +19,4 @@ dependencies:
   - python-dotenv
   - cdsapi>=0.7.6
   - progressbar
+  - scikit-learn
diff --git a/pyproject.toml b/pyproject.toml
index 6229187d..b3f2fb00 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -17,11 +17,12 @@ dependencies = [
     "xarray>=0.20.1",
     "netCDF4",
     "meshio",
-    "geo-skeletons>=0.21.1",
+    "geo-skeletons==0.21.1",
     "geo-parameters>=0.11.0",
     "python-dotenv",
     "progressbar",
     "pyarrow",
+    "scikit-learn",
     "dnplot>=0.4.3",
     "dask[dataframe]"
 ]
diff --git a/tests/test_point_picker/test_point_pickers.py b/tests/test_point_picker/test_point_pickers.py
index 5dd99e37..c367172d 100644
--- a/tests/test_point_picker/test_point_pickers.py
+++ b/tests/test_point_picker/test_point_pickers.py
@@ -1,7 +1,7 @@
 from dnora.grid import Grid
 
 from dnora import pick
-from geo_skeletons import PointSkeleton
+from geo_skeletons import PointSkeleton, GriddedSkeleton
 import numpy as np
 import pytest
 
@@ -58,3 +58,26 @@ def test_selected_points(grid, all_points):
     lat0, lat1 = selected_points.edges("lat")
     assert min(all_points.lat(inds=list(inds))) == lat0
     assert max(all_points.lat(inds=list(inds))) == lat1
+
+
+def test_area_180_wrap():
+    all_points = PointSkeleton(lon=(-175, 179), lat=(20,20))
+    grid = GriddedSkeleton(lon=(-179,-160), lat=(19, 21))
+    inds = pick.Area()(grid=grid, all_points=all_points, expansion_factor=1)
+    assert len(inds) == 1
+    assert inds[0] == 0
+    inds = pick.Area()(grid=grid, all_points=all_points, expansion_factor=2)
+    np.testing.assert_array_almost_equal(inds,[0,1])
+
+def test_nearest_180_wrap():
+    # Nearest point wraps around
+    all_points = PointSkeleton(lon=(-179, 175), lat=(20,20))
+    selected_points = PointSkeleton(lon=(175.1,179.9), lat=(20, 20))
+    inds = pick.NearestGridPoint()(grid=None, all_points=all_points, selected_points=selected_points)
+    np.testing.assert_array_almost_equal(inds,[0,1])
+    
+    # Nearest point is not wrapped around
+    all_points = PointSkeleton(lon=(-149, 175), lat=(20,20))
+    selected_points = PointSkeleton(lon=(175.1,179.9), lat=(20, 20))
+    inds = pick.NearestGridPoint()(grid=None, all_points=all_points, selected_points=selected_points)
+    np.testing.assert_array_almost_equal(inds,[1])
\ No newline at end of file