diff --git a/geodataflow/eogeo/productcatalog.py b/geodataflow/eogeo/productcatalog.py
index 88516df..a9fac06 100644
--- a/geodataflow/eogeo/productcatalog.py
+++ b/geodataflow/eogeo/productcatalog.py
@@ -130,5 +130,5 @@ def _normalize_product(product_type: str, product: object) -> object:
         """
         Normalize properties of specified EOProduct from different EO Providers.
         """
-        setattr(product.geometry, 'srid', 4326)
+        product.geometry = product.geometry.with_srid(4326)
         return product
diff --git a/geodataflow/geoext/commonutils.py b/geodataflow/geoext/commonutils.py
index 652370f..9c214a4 100644
--- a/geodataflow/geoext/commonutils.py
+++ b/geodataflow/geoext/commonutils.py
@@ -39,6 +39,7 @@
 import pyproj as pj
 from shapely.geometry.polygon import Polygon
 from shapely.geometry.base import BaseGeometry
+from shapely.geos import lgeos
 from shapely.ops import transform
 
 from geodataflow.core.capabilities import StoreCapabilities
@@ -354,7 +355,7 @@ def create_transform_function(source_crs: Union[int, str, pj.CRS],
         def transform_fn_(geometry: BaseGeometry) -> BaseGeometry:
             if transform_fn:
                 geometry = transform(transform_fn, geometry)
-                setattr(geometry, 'srid', target_srid)
+                geometry = geometry.with_srid(target_srid)
 
             return geometry
 
@@ -365,12 +366,14 @@ def get_srid(obj: Union[int, object, pj.CRS]) -> int:
         """
         Returns the SRID of this Geometry.
         """
-        if hasattr(obj, 'srid'):
-            return getattr(obj, 'srid')
         if isinstance(obj, int):
             return obj
+        if isinstance(obj, BaseGeometry):
+            return lgeos.GEOSGetSRID(obj._geom)
         if isinstance(obj, pj.CRS):
             return obj.to_epsg()
+        if hasattr(obj, 'srid'):
+            return getattr(obj, 'srid')
 
         return 0
 
@@ -381,7 +384,10 @@ def set_srid(geometry, obj: Union[int, object, pj.CRS]):
         """
         srid = GeometryUtils.get_srid(obj)
         if srid:
-            setattr(geometry, 'srid', srid)
+            if isinstance(geometry, BaseGeometry):
+                lgeos.GEOSSetSRID(geometry._geom, srid)
+            else:
+                setattr(geometry, 'srid', srid)
 
         return geometry
 
diff --git a/geodataflow/geoext/dataset.py b/geodataflow/geoext/dataset.py
index 8f0db0c..bd99174 100644
--- a/geodataflow/geoext/dataset.py
+++ b/geodataflow/geoext/dataset.py
@@ -109,7 +109,7 @@ def geometry(self):
         """
         envelope = self.get_envelope()
         geometry = GeometryUtils.create_geometry_from_bbox(*envelope)
-        setattr(geometry, 'srid', self.get_spatial_srid())
+        geometry = geometry.with_srid(self.get_spatial_srid())
         return geometry
 
     @property
@@ -270,7 +270,7 @@ def clamp_geometry(self, geometry, padding_val: int = 0, minimum_ts: int = 0):
             return envelope
 
         geometry = GeometryUtils.create_geometry_from_bbox(*envelope)
-        setattr(geometry, 'srid', info.get('srid'))
+        geometry = geometry.with_srid(info.get('srid'))
         return geometry
 
     def warp(self,
@@ -285,8 +285,7 @@ def warp(self,
 
         if isinstance(output_geom, list):
             output_geom = GeometryUtils.create_geometry_from_bbox(*output_geom)
-            setattr(output_geom,
-                    'srid', GeometryUtils.get_srid(output_crs) if output_crs else info.get('srid'))
+            output_geom = output_geom.with_srid(GeometryUtils.get_srid(output_crs) if output_crs else info.get('srid'))
 
         if output_geom:
             source_crs = GeometryUtils.get_spatial_crs(info.get('srid'))
diff --git a/geodataflow/geoext/featurestore.py b/geodataflow/geoext/featurestore.py
index 8b0bd1f..a911cf7 100644
--- a/geodataflow/geoext/featurestore.py
+++ b/geodataflow/geoext/featurestore.py
@@ -170,7 +170,7 @@ def features(self) -> Iterable:
             for feature in feature_layer.features():
                 geometry = feature.GetGeometryRef()
                 geometry = shapely_wkb_loads(bytes(geometry.ExportToWkb()))
-                setattr(geometry, 'srid', schema_def.srid)
+                geometry = geometry.with_srid(schema_def.srid)
 
                 feature = type('Feature', (object,), {
                     'type': 'Feature',
diff --git a/geodataflow/pipeline/filters/RasterClip.py b/geodataflow/pipeline/filters/RasterClip.py
index d3cce4c..4a5143f 100644
--- a/geodataflow/pipeline/filters/RasterClip.py
+++ b/geodataflow/pipeline/filters/RasterClip.py
@@ -98,7 +98,7 @@ def run(self, data_store, processing_args):
                 clipping_geoms = [transform_fn(g) for g in clipping_geoms]
             else:
                 for g in clipping_geoms:
-                    setattr(g, 'srid', schema_def.srid)
+                    g = g.with_srid(schema_def.srid)
 
             for dataset in data_store:
                 for g in clipping_geoms:
diff --git a/geodataflow/pipeline/filters/TableUnpack.py b/geodataflow/pipeline/filters/TableUnpack.py
index b8f0c40..990c6b8 100644
--- a/geodataflow/pipeline/filters/TableUnpack.py
+++ b/geodataflow/pipeline/filters/TableUnpack.py
@@ -69,7 +69,7 @@ def run(self, feature_store, processing_args):
             #
             for index, row in data_frame.iterrows():
                 geometry = row['geometry']
-                setattr(geometry, 'srid', schema_def.srid)
+                geometry = geometry.with_srid(schema_def.srid)
 
                 feature = type('Feature', (object,), {
                     'type': 'Feature',