option to reproject data is added #11

openearth · Nov 3, 2020 · 9cf3092 · 9cf3092
1 parent 12a4c3c
commit 9cf3092
Show file tree

Hide file tree

Showing 5 changed files with 25 additions and 8 deletions.
diff --git a/data/osm/background_projected.pgw b/data/osm/background_projected.pgw
@@ -0,0 +1,6 @@
+3.15349231625726656
+0
+0
+-3.15349231625726656
+524565.89883953321259469
+9249556.25685384124517441
diff --git a/data/osm/background_projected.png b/data/osm/background_projected.png
diff --git a/delft3dfmpy/datamodels/common.py b/delft3dfmpy/datamodels/common.py
@@ -87,7 +87,7 @@ def delete_all(self):
             self.iloc[:, 0] = np.nan
             self.dropna(inplace=True)
 
-    def read_shp(self, path, index_col=None, column_mapping=None, check_columns=True, clip=None, check_geotype=True,
+    def read_shp(self, path, index_col=None, column_mapping=None, check_columns=True, proj_crs=None, clip=None, check_geotype=True,
                  id_col='code',filter_cols = False, draintype_col=None, filter_culverts=False, logger=logging):
         """
         Import function, extended with type checks. Does not destroy reference to object.
@@ -115,7 +115,7 @@ def read_shp(self, path, index_col=None, column_mapping=None, check_columns=True
         gdf.drop(gdf.index[gdf.geometry.isnull()], inplace =True) # temporary fix
         logger.debug(f'{missing_features} out of {total_features} do not have a geometry')
 
-        #FIXME: add reprojection of geodataframe via crs inifile
+
 
         if 'MultiPolygon' or 'MultiLineString' in str(gdf.geometry.type):
             #gdf = gdf[gdf.geometry.type != 'MultiPolygon']
@@ -144,6 +144,12 @@ def read_shp(self, path, index_col=None, column_mapping=None, check_columns=True
         if clip is not None:
             self.clip(clip)
 
+        # To re-project CRS system to projected CRS, first all empty geometries should be dropped
+        if proj_crs is not None:
+            gdf.to_crs(epsg=proj_crs, inplace=True)
+        else:
+            logger.info(f'Check if crs of OSM data is an projected crs')
+
 
     def set_data(self, gdf, index_col=None, check_columns=True, check_geotype=True):
 

diff --git a/examples/test_osm.py b/examples/test_osm.py
@@ -35,14 +35,17 @@
 # Id column
 id = config.get('datacolumns','idcolumn')
 
+# Projected coordinate system
+crs_proj = config.get('parameter','projectedcrs')
+
 # TODO: BRANCHES - read id column from json.  Do not deviate between drain type
 # Read branches and store in OSM data model
-osm.branches.read_shp(os.path.join(path,config.get('input','datafile')),index_col=id, clip = osm.clipgeo
+osm.branches.read_shp(os.path.join(path,config.get('input','datafile')),index_col=id, proj_crs= crs_proj, clip = osm.clipgeo
                       , id_col=id, filter_cols=True, logger=logger)
 
 # TODO: CROSS SECTIONS DEFINTION - read id, drain_type, material, width, depth, top_width, diameter, profile_op, profile_cl, bottom_width columns from json
 # read cross-sections
-osm.parametrised_profiles.read_shp(os.path.join(path,config.get('input','datafile')),index_col=id, clip = osm.clipgeo
+osm.parametrised_profiles.read_shp(os.path.join(path,config.get('input','datafile')),index_col=id, proj_crs= crs_proj, clip = osm.clipgeo
                       , id_col=id, filter_cols=True, logger=logger)
 
 # TODO: CROSS SECTIONS DEFINTION - specify roughness dependent on material add this
@@ -65,15 +68,17 @@
 ax.xaxis.set_visible(False)
 ax.yaxis.set_visible(False)
 
-background = plt.imread(path+'/background.png')
-ax.imshow(background, extent=(39.222335471, 39.266427395, -6.814710925, -6.789116518), interpolation='lanczos')
+#background = plt.imread(path+'/background.png')
+#ax.imshow(background, extent=(39.222335471, 39.266427395, -6.814710925, -6.789116518), interpolation='lanczos')
+background = plt.imread(path+'/background_projected.png')
+ax.imshow(background, extent=(524564.3221, 529442.7747, 9246725.9975, 9249557.8336), interpolation='lanczos')
 osm.branches.plot(ax=ax, label='Channel')
 osm.parametrised_profiles.plot(ax=ax, color='C3', label='Cross section')
 plt.show()
 
 # TODO: STRUCTURE - read id, draintype
 # Read culverts
-osm.culverts.read_shp(os.path.join(path,config.get('input','datafile')),index_col=id, clip = osm.clipgeo,
+osm.culverts.read_shp(os.path.join(path,config.get('input','datafile')),index_col=id, proj_crs= crs_proj, clip = osm.clipgeo,
                       id_col=id, filter_cols=True, draintype_col=config.get('datacolumns','draintypecolumn')
                       , filter_culverts=True, logger=logger)
 

diff --git a/notebooks/Usage_introduction_FM_only.ipynb b/notebooks/Usage_introduction_FM_only.ipynb
@@ -764,7 +764,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.8"
+   "version": "3.6.9"
   }
  },
  "nbformat": 4,