Added syst coordinates verification function

modules/hrudelin_1_init.py

@@ -38,11 +38,15 @@ except Exception as e:
 from grass.script.utils import decode, encode
 
 from osgeo import ogr
+from osgeo import gdal
+from osgeo import osr
+
 
 from utils import cutting_raster, get_raster_bounds, get_polygon_bounds 
 from utils import isfloat, isint, get_coords, is_valid_shp, is_valid_geometry
 from utils import is_column_exist, is_column_valid, is_column_value_unique
 from utils import are_columns_value_unique
+from utils import get_Raster_EPSG, get_Vector_EPSG
 
 MY_ABS_PATH=os.path.abspath(__file__)
 MY_DIR=os.path.dirname(MY_ABS_PATH)
@@ -101,20 +105,35 @@ def main(parms_file):
     if not os.path.isfile(dem):
         msg="------------> ERROR : Input Dem "+ dem+" not found"
         sys.exit(msg)
-        
+
+    ## get the ESPG of the DEM
+    demEPSG=get_Raster_EPSG(dem)
+    
+    
     ## if hgeon, landuse and soil provided, test if exist
     for data in parms.items('data'):
         data = os.path.join(directory, data[1])
         if not os.path.isfile(data):
             print(data)
             sys.exit('------------> ERROR : Input data not found' )  
-
+        ## test if EPSG of data is the same of the DEM
+        dataEPSG=get_Raster_EPSG(data)
+        
+        if (demEPSG != dataEPSG):
+            print('input data %s has EPSG %s, input dem %s has EPSG %s'%(data,dataEPSG,dem,demEPSG))
+            sys.exit('------------> ERROR : Input data Syst Coordonnate is not the same than the input dem')  
         
      
     ## Test if gauges exist
     gauges = os.path.join(directory, str(parms.get('gauges', 'gauges')))
     is_valid_shp(gauges)
-       
+
+    ## test if EPSG of gauges is the same of the DEM
+    gaugesEPSG=get_Vector_EPSG(gauges)
+    if (demEPSG != gaugesEPSG):
+            print('input gauges %s has EPSG %s, input dem %s has EPSG %s'%(gauges,gaugesEPSG,dem,demEPSG))
+            sys.exit('------------> ERROR : Input gauges Syst Coordonnate is not the same than the input dem')     
+
     ## Test if gauges is point or multipoint
     is_valid_geometry(gauges,'Point')
             
@@ -143,6 +162,10 @@ def main(parms_file):
     gauges_file = parms.get('gauges', 'relocated_gauges')
     if gauges_file != '':
         is_valid_shp(gauges_file)
+        gauges_fileEPSG=get_Vector_EPSG(gauges_file)
+        if (demEPSG != gauges_fileEPSG):
+            print('input relocated gauges %s has EPSG %s, input dem %s has EPSG %s'%(gauges_file,gauges_fileEPSG,dem,demEPSG))
+            sys.exit('------------> ERROR : Input relocated gauges Syst Coordonnate is not the same than the input dem')     
             
 
          
@@ -154,6 +177,12 @@ def main(parms_file):
         
         ## Test if irrigation is point or multipoint
         is_valid_geometry(irrig,'Point')
+
+        ## test if EPSG of irrigation is the same of the DEM
+        irrigationEPSG=get_Vector_EPSG(irrig)
+        if (demEPSG != irrigationEPSG):
+            print('input irrigation %s has EPSG %s, input dem %s has EPSG %s'%(irrig,irrigationEPSG,dem,demEPSG))
+            sys.exit('------------> ERROR : Input irrigation Syst Coordonnate is not the same than the input dem')    
                 
         ## Test irrig_col_name
         irrig_col_name = parms.get('irrigation', 'irrig_col_name')
@@ -201,6 +230,13 @@ def main(parms_file):
               
         ## Test if irrig_sector is polyon or multipolygon
         is_valid_geometry(irrig_sector,'Polygon')
+
+
+        ## test if EPSG of irrig_sector is the same of the DEM
+        irrig_sectorEPSG=get_Vector_EPSG(irrig_sector)
+        if (demEPSG != irrig_sectorEPSG):
+            print('input irrigation sector %s has EPSG %s, input dem %s has EPSG %s'%(irrig_sector,irrig_sectorEPSG,dem,demEPSG))
+            sys.exit('------------> ERROR : Input irrigation sector Syst Coordonnate is not the same than the input dem')    
             
         ## Test irrig_sector_col_name
         irrig_sector_col_name = parms.get('irrigation', 'irrig_sector_col_name')
@@ -232,6 +268,12 @@ def main(parms_file):
         irrigation_file = parms.get('irrigation', 'relocated_irrigation')
         if irrigation_file != '':
             is_valid_shp(irrigation_file) 
+
+            ## test if EPSG of relocated irrigation is the same of the DEM
+            irrigation_fileEPSG=get_Vector_EPSG(irrigation_file)
+            if (demEPSG != irrigation_fileEPSG):
+                print('input relocated irrigation %s has EPSG %s, input dem %s has EPSG %s'%(irrigation_file,irrigation_fileEPSG,dem,demEPSG))
+                sys.exit('------------> ERROR : Input relocated irrigation Syst Coordonnate is not the same than the input dem')    
             
 
     ## if dams is yes, test input data
@@ -242,7 +284,13 @@ def main(parms_file):
         ## Test if dams is point or multipoint
         is_valid_geometry(dams,'Point')
             
-        
+        ## test if EPSG of dams is the same of the DEM
+        damsEPSG=get_Vector_EPSG(dams)
+        if (demEPSG != damsEPSG):
+            print('input dams %s has EPSG %s, input dem %s has EPSG %s'%(dams,damsEPSG,dem,demEPSG))
+            sys.exit('------------> ERROR : Input dams Syst Coordonnate is not the same than the input dem')    
+
+
         ## Test dams_col_name
         dams_col_name = parms.get('dams', 'dams_col_name')
         if dams_col_name == '':
@@ -283,7 +331,13 @@ def main(parms_file):
 
         dams_file = parms.get('dams', 'relocated_dams')
         if dams_file != '':
-            is_valid_shp(dams_file) 
+            is_valid_shp(dams_file)
+
+            ## test if EPSG of relocated dams is the same of the DEM
+            dams_fileEPSG=get_Vector_EPSG(dams_file)
+            if (demEPSG != dams_fileEPSG):
+                print('input relocated dams %s has EPSG %s, input dem %s has EPSG %s'%(dams_file,dams_fileEPSG,dem,demEPSG))
+                sys.exit('------------> ERROR : Input relocated dams Syst Coordonnate is not the same than the input dem')     
             
 
         #test if dams_col_name and gauges_col_name have no common value

modules/utils.py

@@ -21,6 +21,7 @@ import os, sys
 import pandas as pd
 from osgeo import gdal
 from osgeo import ogr
+from osgeo import osr
 
 def myJoin(f1Path, f2Path, resultPath):
     df1 = pd.read_table(f1Path, delim_whitespace=True)
@@ -180,3 +181,21 @@ def write_log(logf,msg,  start, end, best_x, best_y,point_code,point_type,point_
             'vertic shift': str(abs(vertic_shift))
         })
 
+def get_Raster_EPSG(rasterName):
+    RasterFile=gdal.Open(rasterName)
+    crs=osr.SpatialReference(wkt=RasterFile.GetProjection())
+    if crs:
+        return crs.GetAttrValue('AUTHORITY',1)
+    else:
+        return 0    
+    
+
+def get_Vector_EPSG(vectorName):
+    NameDs=ogr.Open(vectorName)
+    NameLayer=NameDs.GetLayer()
+    crs=NameLayer.GetSpatialRef()
+    if crs:
+        return crs.GetAttrValue('AUTHORITY',1)
+    else:
+        return 0    
+