diff --git a/exposurelib/database.py b/exposurelib/database.py
index 5187211b68a05f2e8e261170716e8dc6d3916981..a58820a39883b3bb2235c5a314203645d8e3732c 100644
--- a/exposurelib/database.py
+++ b/exposurelib/database.py
@@ -875,10 +875,15 @@ class SpatiaLiteExposure(SpatiaLiteDatabase, AbstractExposure):
         self.connection.execute(sql_statement)
         logger.debug("Geometry entry for view %s created" % view_name)
 
-    def import_from_postgis(self, postgis_config, bounding_box=None, country_iso_code=None):
+    def import_from_postgis(
+        self, postgis_config, bounding_box=None, country_iso_code=None, copy_reference=False
+    ):
         """
         Imports exposure data from a PostGIS database limited to either the area of a bounding
-        box or to a country specified by its ISO 3166-1 alpha-3 code.
+        box or to a country specified by its ISO 3166-1 alpha-3 code into the `Asset`, `Entity`
+        and `Taxonomy` tables.
+        If `copy_reference` is set, the function imports reference exposure data (the
+        exposure model created after execution of the Initializer)
         This function first imports all entities while preserving their IDs. Then all assets and
         taxonomies are imported, again preserving the IDs to ensure that the joins between the
         tables do not break.
@@ -896,6 +901,8 @@ class SpatiaLiteExposure(SpatiaLiteDatabase, AbstractExposure):
                 `lon_max`, `lat_min`, `lat_max`.
             country_iso_code (str):
                 ISO 3166-1 alpha-3 code of the country.
+            copy_reference (bool):
+                If set to `True`, the reference exposure model will be imported.
         """
 
         # Empty the tables to be filled
@@ -911,6 +918,16 @@ class SpatiaLiteExposure(SpatiaLiteDatabase, AbstractExposure):
         postgis_db = PostGISDatabase(**postgis_config)
         postgis_db.connect()
 
+        osm_id = "NULL"
+        if copy_reference:
+            input_entity_table = "EntityReference"
+            input_asset_table = "AssetReference"
+            osm_id_extraction_statement = ""
+        else:
+            input_entity_table = "Entity"
+            input_asset_table = "Asset"
+            osm_id_extraction_statement = "Entity.osm_id,"
+
         area_restriction_clause = ""
         # Make sure there is exactly one exposure restricting input
         if (bounding_box and country_iso_code) or (not bounding_box and not country_iso_code):
@@ -918,51 +935,55 @@ class SpatiaLiteExposure(SpatiaLiteDatabase, AbstractExposure):
 
         # Create the exposure import-area restrictive SQL statement
         if bounding_box is not None:
-            area_restriction_clause = """ WHERE geom && ST_MakeEnvelope(%f, %f, %f, %f, 4326)
-                """ % (
-                bounding_box["lon_min"],
-                bounding_box["lat_min"],
-                bounding_box["lon_max"],
-                bounding_box["lat_max"],
-            )
-        elif country_iso_code is not None:
-            area_restriction_clause = (
+            area_restriction_clause = f""" WHERE geom && ST_MakeEnvelope
+                (
+                {bounding_box["lon_min"]},
+                {bounding_box["lat_min"]},
+                {bounding_box["lon_max"]},
+                {bounding_box["lat_max"]},
+                4326
+                )
                 """
-                INNER JOIN obm_tiles ON obm_tiles.quadkey = Entity.quadkey
-                WHERE obm_tiles.country_iso_code = '%s'
+        elif country_iso_code is not None:
+            area_restriction_clause = f"""
+                INNER JOIN obm_tiles ON obm_tiles.quadkey = {input_entity_table}.quadkey
+                WHERE obm_tiles.country_iso_code = '{country_iso_code}'
                 """
-                % country_iso_code
-            )
 
         # Import entities
-        sql_statement = (
+        sql_statement = f"""
+            SELECT {input_entity_table}.id, {input_entity_table}.quadkey,
+            {osm_id_extraction_statement} ST_AsText({input_entity_table}.geom)
+            FROM {input_entity_table}
+            {area_restriction_clause}
             """
-            SELECT Entity.id, Entity.quadkey, Entity.osm_id, ST_AsText(Entity.geom) FROM Entity
-            %s """
-            % area_restriction_clause
-        )
+
         postgis_db.cursor.execute(sql_statement)
         entities = postgis_db.cursor.fetchall()
-        for entity_id, quadkey, osm_id, wkt_geom in entities:
+        for entity_info in entities:
+            if copy_reference:
+                entity_id, quadkey, wkt_geom = entity_info
+            else:
+                entity_id, quadkey, osm_id, wkt_geom = entity_info
             if not osm_id:
                 osm_id = "NULL"
             sql_statement = f"""
                 INSERT INTO Entity
                 (id, quadkey, osm_id, geom)
                 VALUES ({entity_id}, '{quadkey}', {osm_id},
-                    ST_Multi(ST_GeomFromText('{wkt_geom}', 4326)))
+                ST_Multi(ST_GeomFromText('{wkt_geom}', 4326)))
                 """
             self.cursor.execute(sql_statement)
 
         # Import assets
-        sql_statement = (
-            """
+        sql_statement = f"""
             SELECT entity_id, taxonomy_id, number, structural, night
-            FROM Asset
-            INNER JOIN Entity
-            ON Entity.id = Asset.entity_id %s """
-            % area_restriction_clause
-        )
+            FROM {input_asset_table}
+            INNER JOIN {input_entity_table}
+            ON
+            {input_entity_table}.id = {input_asset_table}.entity_id
+            {area_restriction_clause}
+            """
         postgis_db.cursor.execute(sql_statement)
         assets = postgis_db.cursor.fetchall()
         for entity_id, taxonomy_id, number, structural, night in assets: