diff --git a/climada/util/coordinates.py b/climada/util/coordinates.py
index 3b49c74b8..23d1ce0d9 100644
--- a/climada/util/coordinates.py
+++ b/climada/util/coordinates.py
@@ -2412,13 +2412,16 @@ def apply_box(df_exp):
     LOGGER.info('Raster from resolution %s to %s.', res, raster_res)
     df_poly = gpd.GeoDataFrame(points_df[val_names])
     if not scheduler:
-        df_poly['geometry'] = apply_box(points_df)
+        df_poly['_-geometry-prov'] = apply_box(points_df)
     else:
         ddata = dd.from_pandas(points_df[['latitude', 'longitude']],
                                npartitions=cpu_count())
-        df_poly['geometry'] = ddata.map_partitions(apply_box, meta=Polygon) \
+        df_poly['_-geometry-prov'] = ddata.map_partitions(apply_box, meta=Polygon) \
                                    .compute(scheduler=scheduler)
-    df_poly.set_crs(crs if crs else points_df.crs if points_df.crs else DEF_CRS, inplace=True)
+    df_poly.set_geometry('_-geometry-prov',
+                         crs=crs if crs else points_df.crs if points_df.crs else DEF_CRS,
+                         inplace=True,
+                         drop=True)
 
     # renormalize longitude if necessary
     if equal_crs(df_poly.crs, DEF_CRS):
@@ -2651,11 +2654,7 @@ def set_df_geometry_points(df_val, scheduler=None, crs=None):
     LOGGER.info('Setting geometry points.')
 
     # keep the original crs if any
-    if crs is None:
-        try:
-            crs = df_val.geometry.crs
-        except AttributeError:
-            crs = None
+    crs = df_val.crs if crs is None else crs  # crs might now still be None
 
     # work in parallel
     if scheduler:
@@ -2663,18 +2662,20 @@ def apply_point(df_exp):
             return df_exp.apply(lambda row: Point(row.longitude, row.latitude), axis=1)
 
         ddata = dd.from_pandas(df_val, npartitions=cpu_count())
-        df_val['geometry'] = ddata.map_partitions(
+        df_val['_-geometry-prov'] = ddata.map_partitions(
                                  apply_point,
                                  meta=('geometry', gpd.array.GeometryDtype)
                              ).compute(scheduler=scheduler)
+
     # single process
     else:
-        df_val['geometry'] = gpd.GeoSeries(
-            gpd.points_from_xy(df_val.longitude, df_val.latitude), index=df_val.index, crs=crs)
+        df_val['_-geometry-prov'] = gpd.GeoSeries(
+            gpd.points_from_xy(df_val.longitude, df_val.latitude),
+            index=df_val.index)
 
-    # set crs
-    if crs:
-        df_val.set_crs(crs, inplace=True)
+    # A 'geometry' column must not be created in a GeoDataFrame except through the constructor
+    # or with set_geometry. That's why we first made a temporary columns with a weird name
+    df_val.set_geometry('_-geometry-prov', inplace=True, drop=True, crs=crs)
 
 
 def fao_code_def():