Crs extraction

.travis.yml

@@ -3,7 +3,6 @@ dist: xenial
 language: python
 python:
   - "3.8"
-  - "nightly"
 addons:
   apt:
     sources:

geoextent/__init__.py

@@ -1,3 +1,3 @@
 name = "geoextent"
 
-__version__ = '0.6.0'
+__version__ = '0.6.1'

geoextent/lib/extent.py

@@ -34,6 +34,18 @@ def compute_bbox_wgs84(module, path):
         raise Exception("The bounding box could not be transformed to the target CRS epsg:{} \n error {}"
                         .format(hf.WGS84_EPSG_ID, e))
 
+    validate = hf.validate_bbox_wgs84(spatial_extent['bbox'])
+    logger.debug("Validate: {}".format(validate))
+
+    if not hf.validate_bbox_wgs84(spatial_extent['bbox']):
+        try:
+            flip_bbox = hf.flip_bbox(spatial_extent['bbox'])
+            spatial_extent['bbox'] = flip_bbox
+
+        except Exception as e:
+            raise Exception("The bounding box could not be transformed to the target CRS epsg:{} \n error {}"
+                            .format(hf.WGS84_EPSG_ID, e))
+
     return spatial_extent
 
 

geoextent/lib/handleCSV.py

@@ -3,7 +3,6 @@
 from osgeo import gdal
 from . import helpfunctions as hf
 
-
 logger = logging.getLogger("geoextent")
 
 search = {"longitude": ["(.)*longitude", "(.)*long(.)*", "^lon", "lon$", "(.)*lng(.)*", "^x", "x$"],
@@ -30,8 +29,8 @@ def checkFileSupported(filepath):
 
     if driver == "CSV":
         with open(filepath) as csv_file:
-            daten = csv.reader(csv_file.readlines())
-            if daten is None:
+            data = csv.reader(csv_file.readlines())
+            if data is None:
                 logger.debug("File {} is NOT supported by HandleCSV module".format(filepath))
                 return False
             else:
@@ -40,6 +39,7 @@ def checkFileSupported(filepath):
     else:
         return False
 
+
 def getBoundingBox(filePath):
     '''
     Function purpose: extracts the spatial extent (bounding box) from a csv-file \n
@@ -49,32 +49,32 @@ def getBoundingBox(filePath):
 
     with open(filePath) as csv_file:
         # To get delimiter either comma or simecolon
-        daten = hf.getDelimiter(csv_file)
+        date = hf.getDelimiter(csv_file)
 
         elements = []
-        for x in daten:
+        for x in date:
             elements.append(x)
 
-        spatialLatExtent = hf.searchForParameters(elements, search['latitude'], exp_data='numeric')
+        spatial_lat_extent = hf.searchForParameters(elements, search['latitude'], exp_data='numeric')
 
-        minlat = None
-        maxlat = None
+        min_lat = None
+        max_lat = None
 
-        if spatialLatExtent is None:
+        if spatial_lat_extent is None:
             pass
         else:
-            minlat = (min(spatialLatExtent))
-            maxlat = (max(spatialLatExtent))
+            min_lat = (min(spatial_lat_extent))
+            max_lat = (max(spatial_lat_extent))
 
-        spatialLonExtent = hf.searchForParameters(elements, search['longitude'], exp_data='numeric')
+        spatial_lon_extent = hf.searchForParameters(elements, search['longitude'], exp_data='numeric')
 
-        if spatialLonExtent is None:
+        if spatial_lon_extent is None:
             raise Exception('The csv file from ' + filePath + ' has no BoundingBox')
         else:
-            minlon = (min(spatialLonExtent))
-            maxlon = (max(spatialLonExtent))
+            min_lon = (min(spatial_lon_extent))
+            max_lon = (max(spatial_lon_extent))
 
-        bbox = [minlon, minlat, maxlon, maxlat]
+        bbox = [min_lon, min_lat, max_lon, max_lat]
         logger.debug("Extracted Bounding box (without projection): {}".format(bbox))
         crs = getCRS(filePath)
         logger.debug("Extracted CRS: {}".format(crs))
@@ -84,53 +84,59 @@ def getBoundingBox(filePath):
 
     return spatialExtent
 
-def getTemporalExtent(filePath, num_sample):
-    ''' extract time extent from csv string \n
+
+def getTemporalExtent(filepath, num_sample):
+    """ extract time extent from csv string \n
     input "filePath": type string, file path to csv File \n
     returns temporal extent of the file: type list, length = 2, both entries have the type str, temporalExtent[0] <= temporalExtent[1]
-    '''
+    """
 
-    with open(filePath) as csv_file:
-        # To get delimiter either comma or simecolon
-        daten = hf.getDelimiter(csv_file)
+    with open(filepath) as csv_file:
+        # To get delimiter either comma or semicolon
+        data = hf.getDelimiter(csv_file)
 
         elements = []
-        for x in daten:
+        for x in data:
             elements.append(x)
 
         all_temporal_extent = hf.searchForParameters(elements, search['time'], exp_data="time")
         if all_temporal_extent is None:
-            raise Exception('The csv file from ' + filePath + ' has no TemporalExtent')
+            raise Exception('The csv file from ' + filepath + ' has no TemporalExtent')
         else:
             tbox = []
-            parsed_time = hf.date_parser(all_temporal_extent, num_sample = num_sample)
+            parsed_time = hf.date_parser(all_temporal_extent, num_sample=num_sample)
 
             if parsed_time is not None:
                 # Min and max into ISO8601 format ('%Y-%m-%d')
                 tbox.append(min(parsed_time).strftime('%Y-%m-%d'))
                 tbox.append(max(parsed_time).strftime('%Y-%m-%d'))
             else:
-                raise Exception('The csv file from ' + filePath + ' has no recognizable TemporalExtent')
+                raise Exception('The csv file from ' + filepath + ' has no recognizable TemporalExtent')
             return tbox
 
-def getCRS(filePath):
+
+def getCRS(filepath):
     '''extracts coordinatesystem from csv File \n
     input "filepath": type string, file path to csv file \n
     returns the epsg code of the used coordinate reference system, type list, contains extracted coordinate system of content from csv file
     '''
 
-    with open(filePath) as csv_file:
-        daten = csv.reader(csv_file.readlines())
+    with open(filepath) as csv_file:
+        data = hf.getDelimiter(csv_file)
         elements = []
-        for x in daten:
+        for x in data:
             elements.append(x)
 
-        if hf.searchForParameters(elements, search['latitude'] + search['longitude']) is None:
-            if hf.searchForParameters(elements, ["crs","srsID"]) is None:
-                raise Exception('The csv file from ' + filePath + ' has no CRS')
-            if hf.searchForParameters(elements, ["crs","srsID"]) == "WGS84":
-                return "4326"
-            else:
-                raise Exception('The csv file from ' + filePath + ' has no WGS84 CRS')
+        param = hf.searchForParameters(elements, ["crs", "srsID", "EPSG"])
+
+        if param is None:
+            logger.debug("{} : There is no identifiable coordinate reference system. We will try to use EPSG: 4326".format(filepath))
+            crs = "4326"
+
+        elif len(list(set(param))) > 1:
+            logger.debug("{} : Coordinate reference system of the file is ambiguous. Extraction is not possible.".format(filepath))
+            raise Exception('The csv file from ' + filepath + ' has no CRS')
         else:
-            return "4326"
+            crs = str(list(set(param))[0])
+
+        return crs

geoextent/lib/handleRaster.py

@@ -33,17 +33,17 @@ def checkFileSupported(filepath):
         return False
 
 
-def getBoundingBox(filePath):
-    ''' extracts bounding box from geotiff \n
-    input "filepath": type string, file path to geotiff file \n
-    returns bounding box of the file: type list, length = 4 , type = float, schema = [min(longs), min(lats), max(longs), max(lats)] 
-    '''
+def getBoundingBox(filepath):
+    """ extracts bounding box from raster \n
+    input "filepath": type string, file path to raster file \n
+    returns bounding box of the file: type list, length = 4 , type = float, schema = [min(longs), min(lats), max(longs), max(lats)]
+    """
     # Enable exceptions
 
     crs_output = hf.WGS84_EPSG_ID
     gdal.UseExceptions()
 
-    geotiffContent = gdal.Open(filePath)
+    geotiffContent = gdal.Open(filepath)
 
     # get the existing coordinate system
     old_crs = osr.SpatialReference()
@@ -59,33 +59,38 @@ def getBoundingBox(filePath):
     height = geotiffContent.RasterYSize
     gt = geotiffContent.GetGeoTransform()
 
-    minx = gt[0]
-    miny = gt[3] + width * gt[4] + height * gt[5]
-    maxx = gt[0] + width * gt[1] + height * gt[2]
-    maxy = gt[3]
+    min_x = gt[0]
+    min_y = gt[3] + width * gt[4] + height * gt[5]
+    max_x = gt[0] + width * gt[1] + height * gt[2]
+    max_y = gt[3]
 
     transform = osr.CoordinateTransformation(old_crs, new_crs)
-    # get the coordinates in lat long
-    latlongmin = transform.TransformPoint(minx, miny)
-    latlongmax = transform.TransformPoint(maxx, maxy)
+    try:
+        # get the coordinates in lat long
+        lat_long_min = transform.TransformPoint(min_x, min_y)
+        lat_long_max = transform.TransformPoint(max_x, max_y)
+    except:
+        # Assume that coordinates are in EPSG:4236
+        logger.debug("{}: There is no identifiable coordinate reference system. We will try to use EPSG: 4326 "
+                     .format(filepath))
+        lat_long_min = [min_y, min_x]
+        lat_long_max = [max_y, max_x]
 
-    bbox = [latlongmin[0], latlongmin[1], latlongmax[0], latlongmax[1]]
+    bbox = [lat_long_min[0], lat_long_min[1], lat_long_max[0], lat_long_max[1]]
 
     if int(osgeo.__version__[0]) >= 3:
         if old_crs.GetAxisMappingStrategy() == 1:
-            bbox = [latlongmin[1], latlongmin[0], latlongmax[1], latlongmax[0]]
+            bbox = [lat_long_min[1], lat_long_min[0], lat_long_max[1], lat_long_max[0]]
 
     spatialExtent = {"bbox": bbox, "crs": str(crs_output)}
 
     return spatialExtent
 
 
-def getTemporalExtent(filePath):
-    ''' extracts temporal extent of the geotiff \n
+def getTemporalExtent(filepath):
+    """ extracts temporal extent of the geotiff \n
     input "filepath": type string, file path to geotiff file \n
-    returns None as There is no time value for GeoTIFF files 
-    '''
-
-    print('There is no time value for GeoTIFF files')
+    returns None as There is no time value for GeoTIFF files
+    """
+    logger.debug('{} There is no time value for raster files'.format(filepath))
     return None
-

geoextent/lib/handleVector.py

@@ -110,15 +110,16 @@ def getBoundingBox(filepath):
     """
     datasource = ogr.Open(filepath)
     geo_dict = {}
-    crs_output = hf.WGS84_EPSG_ID
 
     for layer in datasource:
         layer_name = layer.GetDescription()
         ext = layer.GetExtent()
         bbox = [ext[0], ext[2], ext[1], ext[3]]
 
         try:
-            crs = layer.GetSpatialRef().GetAttrValue("GEOGCS|AUTHORITY", 1)
+            spatial_ref = layer.GetSpatialRef()
+            spatial_ref.AutoIdentifyEPSG()
+            crs = spatial_ref.GetAuthorityCode(None)
         except Exception as e:
             logger.debug("Error extracting EPSG CODE from layer {}: \n {}".format(layer_name, e))
             crs = None

geoextent/lib/helpfunctions.py

@@ -87,7 +87,7 @@ def searchForParameters(elements, paramArray, exp_data=None):
 
 def transformingIntoWGS84(crs, coordinate):
     '''
-    Function purpose: transforming SRS into WGS84 (EPSG:4978; used by the GPS satellite navigation system) \n
+    Function purpose: transforming SRS into WGS84 (EPSG:4326) \n
     Input: crs, point \n
     Output: retPoint constisting of x2, y2 (transformed points)
     '''
@@ -112,7 +112,7 @@ def transformingIntoWGS84(crs, coordinate):
 
 def transformingArrayIntoWGS84(crs, pointArray):
     '''
-    Function purpose: transforming SRS into WGS84 (EPSG:4978; used by the GPS satellite navigation system) from an array \n
+    Function purpose: transforming SRS into WGS84 (EPSG 4326; used by the GPS satellite navigation system) from an array \n
     Input: crs, pointArray \n
     Output: array array
     '''
@@ -130,6 +130,40 @@ def transformingArrayIntoWGS84(crs, pointArray):
         return [transf_bbox[0][0], transf_bbox[0][1], transf_bbox[1][0], transf_bbox[1][1]]
 
 
+def validate_bbox_wgs84(bbox):
+    """
+    :param bbox:
+    :return:
+    """
+    valid = True
+    lon_values = bbox[0:3:2]
+    lat_values = bbox[1:4:2]
+
+    if sum(list(map(lambda x: x < -90 or x > 90, lat_values))) + sum(
+            list(map(lambda x: x < -180 or x > 180, lon_values))) > 0:
+        valid = False
+
+    return valid
+
+
+def flip_bbox(bbox):
+    """
+    :param bbox:
+    :return:
+    """
+    # Flip values
+    lon_values = bbox[1:4:2]
+    lat_values = bbox[0:3:2]
+
+    bbox_flip = [lon_values[0], lat_values[0], lon_values[1], lat_values[1]]
+    if validate_bbox_wgs84(bbox_flip):
+        logger.warning("Longitude and latitude values flipped")
+        return bbox_flip
+    else:
+        raise Exception("Latitude and longitude values extracted do not seem to be correctly transformed. We tried "
+                        "flipping latitude and longitude values but both bbox are incorrect")
+
+
 def validate(date_text):
     try:
         if datetime.datetime.strptime(date_text, output_time_format):