Python osgeo.osr.SpatialReference() Examples

def writeSRS(self, srs):
        """
        @summary: Writes spatial reference system information to this raster file.
        @param srs: An osgeo.osr.SpatialReference object or
                        A WKT string describing the desired spatial reference system.  
                          Raster.populateStats will populate the Raster.srs 
                          attribute with a correctly formatted string.  Correctly 
                          formatted strings are also output by:
                             * osgeo.gdal.Dataset.GetProjection
                             * osgeo.osr.SpatialReference.ExportToWKT 
        @postcondition: The raster file is updated with new srs information.
        @raise LMError: on failure to open dataset or write srs
        """
        from LmServer.common.geotools import GeoFileInfo

        if (self._dlocation is not None and os.path.exists(self._dlocation)):
            geoFI = GeoFileInfo(self._dlocation, updateable=True)
            if isinstance(srs, osr.SpatialReference):
                srs = srs.ExportToWkt()
            geoFI.writeWktSRS(srs)
                
# ............................................... 

def _transformBBox(self, origEpsg=None, origBBox=None):
        if origEpsg is None:
            origEpsg = self._epsg
        if origBBox is None:
            origBBox = self._bbox
        minX, minY, maxX, maxY = origBBox
            
        if origEpsg != DEFAULT_EPSG:
            srcSRS = osr.SpatialReference()
            srcSRS.ImportFromEPSG(origEpsg)
            dstSRS = osr.SpatialReference()
            dstSRS.ImportFromEPSG(DEFAULT_EPSG)
            
            spTransform = osr.CoordinateTransformation(srcSRS, dstSRS)
            # Allow for return of either (x, y) or (x, y, z) 
            retvals = spTransform.TransformPoint(minX, minY)
            newMinX, newMinY = retvals[0], retvals[1]
            retvals = spTransform.TransformPoint(maxX, maxY)
            newMaxX, newMaxY = retvals[0], retvals[1]
            return (newMinX, newMinY, newMaxX, newMaxY)
        else:
            return origBBox
        
# ............................................... 

def __init__(self, lon, lat):
    '''
    initialize projection dataset based on (regular) lat/lon vectors
    '''    
    epsg = 4326 # EPSG code for regular lat/long grid
    # size of dataset
    size = (len(lon), len(lat))
    # GDAL geotransform vector
    dx = lon[1]-lon[0]; dy = lat[1]-lat[0]
    ulx = lon[0]-dx/2.; uly = lat[0]-dy/2. # coordinates of upper left corner (same for source and sink)
    # GT(2) & GT(4) are zero for North-up; GT(1) & GT(5) are pixel width and height; (GT(0),GT(3)) is the top left corner
    geotransform = (ulx, dx, 0., uly, 0., dy) 
    # GDAL projection 
    projection = osr.SpatialReference()
    projection.ImportFromEPSG(epsg)
    # create GDAL projection object from parent instance
    super(LatLonProj,self).__init__(projection=projection, geotransform=geotransform, size=size)
    self.epsg = epsg # save projection code number
    
## function to reproject and resample a 2D array 

def __setstate__(self, pickle):
    ''' support pickling, necessary for multiprocessing: GDAL is not pickable '''
    # handle projection
    self.projection = osr.SpatialReference() 
    self.projection.SetWellKnownGeogCS('WGS84')           
    self.projection.ImportFromWkt(pickle['_projection'])  # from Well-Known-Text
    del pickle['_projection'] # not actually an attribute
    # handle axes
    self.geotransform = pickle['_geotransform']
    self.isProjected = pickle['_isProjected']
    xlon, ylat = getAxes(geotransform=self.geotransform, xlen=pickle['_xlon'], ylen=pickle['_ylat'], 
                         projected=self.isProjected)
    self.xlon = xlon; self.ylat = ylat
    del pickle['_geotransform'], pickle['_isProjected'], pickle['_xlon'], pickle['_ylat']
    # update instance dict with pickle dict
    self.__dict__.update(pickle) 

def __init__(self, back_ds, wkt_stored, rect, **kwargs):
        wkt_virtual = back_ds.virtual_of_stored_given_mode(
            wkt_stored, back_ds.wkt_work, back_ds.wkt_fallback, back_ds.wkt_forced,
        )

        if wkt_virtual is not None:
            sr_virtual = osr.SpatialReference(wkt_virtual)
        else:
            sr_virtual = None

        to_work, to_virtual = back_ds.get_transforms(sr_virtual, rect)

        self.back_ds = back_ds
        self.wkt_stored = wkt_stored
        self.wkt_virtual = wkt_virtual
        self.to_work = to_work
        self.to_virtual = to_virtual
        super().__init__(**kwargs) 

def __init__(self, wkt_work, wkt_fallback, wkt_forced, analyse_transformation, **kwargs):

        if wkt_work is not None:
            sr_work = osr.SpatialReference(wkt_work)
        else:
            sr_work = None
        if wkt_fallback is not None:
            sr_fallback = osr.SpatialReference(wkt_fallback)
        else:
            sr_fallback = None
        if wkt_forced is not None:
            sr_forced = osr.SpatialReference(wkt_forced)
        else:
            sr_forced = None

        self.wkt_work = wkt_work
        self.wkt_fallback = wkt_fallback
        self.wkt_forced = wkt_forced
        self.sr_work = sr_work
        self.sr_fallback = sr_fallback
        self.sr_forced = sr_forced
        self.analyse_transformations = analyse_transformation
        super(BackDatasetConversionsMixin, self).__init__(**kwargs) 

def GetGeoInfo(FileName):

    if exists(FileName) is False:
            raise Exception('[Errno 2] No such file or directory: \'' + FileName + '\'')    
    
    
    SourceDS = gdal.Open(FileName, gdal.GA_ReadOnly)
    if SourceDS == None:
        raise Exception("Unable to read the data file")
    
    NDV = SourceDS.GetRasterBand(1).GetNoDataValue()
    xsize = SourceDS.RasterXSize
    ysize = SourceDS.RasterYSize
    GeoT = SourceDS.GetGeoTransform()
    Projection = osr.SpatialReference()
    Projection.ImportFromWkt(SourceDS.GetProjectionRef())
    DataType = SourceDS.GetRasterBand(1).DataType
    DataType = gdal.GetDataTypeName(DataType)
    
    return NDV, xsize, ysize, GeoT, Projection, DataType
#==============================================================================

#==============================================================================
# Function to read the original file's projection: 

def polygonize(img,shp_path):
    # mapping between gdal type and ogr field type
    type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                    gdal.GDT_UInt16: ogr.OFTInteger,
                    gdal.GDT_Int16: ogr.OFTInteger,
                    gdal.GDT_UInt32: ogr.OFTInteger,
                    gdal.GDT_Int32: ogr.OFTInteger,
                    gdal.GDT_Float32: ogr.OFTReal,
                    gdal.GDT_Float64: ogr.OFTReal,
                    gdal.GDT_CInt16: ogr.OFTInteger,
                    gdal.GDT_CInt32: ogr.OFTInteger,
                    gdal.GDT_CFloat32: ogr.OFTReal,
                    gdal.GDT_CFloat64: ogr.OFTReal}

    ds = gdal.Open(img)
    prj = ds.GetProjection()
    srcband = ds.GetRasterBand(1)
    dst_layername = "Shape"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference(wkt=prj)

    dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
    raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# convert images in a selected folder to shapefiles 

def get_coord_transform(source_epsg, target_epsg):
    '''
    Creates an OGR-framework coordinate transformation for use in projecting
    coordinates to a new coordinate reference system (CRS). Used as, e.g.:
        transform = get_coord_transform(source_epsg, target_epsg)
        transform.TransformPoint(x, y)
    Arguments:
        source_epsg     The EPSG code for the source CRS
        target_epsg     The EPSG code for the target CRS
    '''
    # Develop a coordinate transformation, if desired
    transform = None
    source_ref = osr.SpatialReference()
    target_ref = osr.SpatialReference()
    source_ref.ImportFromEPSG(source_epsg)
    target_ref.ImportFromEPSG(target_epsg)
    return osr.CoordinateTransformation(source_ref, target_ref) 

def get_idx_as_shp(self, path, gt, wkt):
        '''
        Exports a Shapefile containing the locations of the extracted
        endmembers. Assumes the coordinates are in decimal degrees.
        '''
        coords = pixel_to_xy(self.get_idx(), gt=gt, wkt=wkt, dd=True)

        driver = ogr.GetDriverByName('ESRI Shapefile')
        ds = driver.CreateDataSource(path)
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(4326)

        layer = ds.CreateLayer(path.split('.')[0], srs, ogr.wkbPoint)
        for pair in coords:
            feature = ogr.Feature(layer.GetLayerDefn())

            # Create the point from the Well Known Text
            point = ogr.CreateGeometryFromWkt('POINT(%f %f)' %  pair)
            feature.SetGeometry(point) # Set the feature geometry
            layer.CreateFeature(feature) # Create the feature in the layer
            feature.Destroy() # Destroy the feature to free resources

        # Destroy the data source to free resources
        ds.Destroy() 

def test_rasterreader_bboxrounding(classraster_filepath):
    reader = RasterReader(classraster_filepath)
    assert reader
    assert isinstance(reader.srs, osr.SpatialReference)
    assert reader.bbox == Bbox(727000.0, 6171000.0, 728000.0, 6172000.0)
    pix_win = reader.bbox_to_pixel_window(reader.bbox)
    assert pix_win == (0, 0, 500, 500)
    # Try bboxes which are slightly off pixel borders
    pix_win = reader.bbox_to_pixel_window(
        Bbox(727001.9999, 6171000, 728000.0, 6172000.0)
    )
    assert pix_win == (1, 0, 499, 500)
    pix_win = reader.bbox_to_pixel_window(Bbox(727000, 6171000, 727999.99, 6171999.99))
    assert pix_win == (0, 0, 500, 500)
    pix_win = reader.bbox_to_pixel_window(
        Bbox(727000, 6171000, 728000.00001, 6172000.00001)
    )
    assert pix_win == (0, 0, 500, 500) 

def get_map_bounds(mapname):
	if mapname in maps:
		thismap = maps[mapname]
	else:
		print("Map", mapname, "does not exist.")
		return
	xsize, ysize = thismap.RasterXSize, thismap.RasterYSize
	gt = thismap.GetGeoTransform()
	proj = osr.SpatialReference()
	proj.ImportFromWkt(thismap.GetProjection())
	transform = osr.CreateCoordinateTransformation(proj, wgs)
	minp = gm.transform(gt, (0, 0))
	maxp = gm.transform(gt, (xsize, ysize))
	xmin, ymin, _ = transform.TransformPoint(minp[0], minp[1])
	xmax, ymax, _ = transform.TransformPoint(maxp[0], maxp[1])
	north = max(ymin, ymax) # Maps might be reversed, so we're not sure which one is actually the maximum
	east = max(xmin, xmax)
	south = min(ymin, ymax)
	west = min(xmin, xmax)
	return (north, east, south, west) 

def to_projection(self, variable, projection):
        """Convert Grid to New Projection.

            Parameters
            ----------
            variable: :obj:`str`
                Name of variable in dataset.
            projection: :func:`osr.SpatialReference`
                Projection to convert data to.

            Returns
            -------
            :func:`xarray.Dataset`
        """
        new_data = []
        for band in range(self._obj.dims[self.time_dim]):
            arr_grid = ArrayGrid(in_array=self._obj[variable][band].values,
                                 wkt_projection=self.projection.ExportToWkt(),
                                 geotransform=self.geotransform)
            ggrid = arr_grid.to_projection(projection, gdalconst.GRA_Average)
            new_data.append(ggrid.np_array())

        self.to_datetime()
        return self._export_dataset(variable, np.array(new_data),
                                    ggrid) 

def getSpatialReference(self, output, datatype):
        """
        :param output: :class:`pywps.Process.InAndOutputs.ComplexOutput`
        :param datatype: String
        :return: projection of the output
        """

        sr = osr.SpatialReference()
        if datatype == "raster":
            wkt = self.dataset.GetProjection()
            res = sr.ImportFromWkt(wkt)
            if res == 0:
                return sr
        elif datatype == "vector":
            layer = self.dataset.GetLayer()
            ref = layer.GetSpatialRef()
            if ref:
                return ref
        return None 

def get_layer_dimensions(template_layer_filename):
    """
    @summary: Get layer information
    @param template_layer_filename: The template layer to get information for
    """
    ds = gdal.Open(template_layer_filename)

    num_cols = ds.RasterXSize
    num_rows = ds.RasterYSize
    
    minx, cell_size, _, maxy, _, _ = ds.GetGeoTransform()
    
    maxx = minx + (cell_size * num_cols)
    miny = maxy - (cell_size * num_rows)
    
    bbox = (minx, miny, maxx, maxy)

    srs = osr.SpatialReference()
    srs.ImportFromWkt(ds.GetProjectionRef())
    # TODO: Try to find a better way to fail
    try:
        epsg = int(srs.GetAttrValue('AUTHORITY', 1))
    except:
        epsg = DEFAULT_EPSG
    ds = None
    
    return bbox, cell_size, epsg

# ............................................................................. 

def write_tiff(out_filename, bbox, cell_size, data, epsg,
               nodata=DEFAULT_NODATA):
    """
    @summary: Write a GeoTiff raster layer from a numpy array
    @param out_filename: The file location to write the raster to
    @param bbox: (minx, miny, maxx, maxy) tuple of raster coordinates
    @param cell_size: The cell size, in map units, of each cell in the raster
    @param data: The data array to use for raster content
    @param epsg: The epsg code of the map projection to use for this raster
    @param nodata: A value to use for NODATA
    @todo: Probably should establish a common layer package in LmCommon
    """
    minx, miny, maxx, maxy = bbox
    (num_rows, num_cols) = data.shape
    
    readyFilename(out_filename, overwrite=True)
    
    drv = gdal.GetDriverByName(LMFormat.GTIFF.driver)
    ds = drv.Create(out_filename, num_cols, num_rows, 1, gdalconst.GDT_Byte)
    ds.SetGeoTransform([minx, cell_size, 0, maxy, 0, -cell_size])

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg)
    ds.SetProjection(srs.ExportToWkt())
    
    outBand = ds.GetRasterBand(1)
    outBand.WriteArray(data)
    outBand.FlushCache()
    outBand.SetNoDataValue(nodata)

# .............................................................................
# Testing code 

def getSRS(self):
        if (self._dlocation is not None and os.path.exists(self._dlocation)):
            ds = gdal.Open(str(self._dlocation), gdalconst.GA_ReadOnly)
            wktSRS = ds.GetProjection()
            if wktSRS is not '':
                srs = osr.SpatialReference()
                srs.ImportFromWkt(wktSRS)
            else:
                srs = self.createSRSFromEPSG()
            ds = None
            return srs
        else:
            raise LMError(currargs='Input file %s does not exist' % self._dlocation)

# ............................................... 

def createSRSFromEPSG(self, epsgcode=None):
        srs = SpatialReference()
        if epsgcode is not None:
            srs.ImportFromEPSG(epsgcode)
        else:
            srs.ImportFromEPSG(self._epsg)
        return srs

# ............................................... 

def translatePoints(points, srcWKT=None, srcEPSG=None, dstWKT=None, dstEPSG=None):
        """
        @summary: Translates a list of (x, y) pairs from one coordinate system to 
                         another
        @param points: A list of points [(x, y)+]
        @param srcWKT: The Well-Known Text for the source coordinate system
        @param srcEPSG: The EPSG code of the source coordinate system
        @param dstWKT: The Well-Known Text for the destination coordinate system
        @param dstEPSG: The EPSG code of the destination coordinate system
        @note: Use either srcWKT or srcEPSG, srcWKT has priority
        @note: Use either dstWKT or dstEPSG, dstWKT has priority
        """
        srcSR = SpatialReference()
        if srcWKT is not None:
            srcSR.ImportFromWkt(srcWKT)
        elif srcEPSG is not None:
            srcSR.ImportFromEPSG(srcEPSG)
        else:
            raise Exception("Either srcWKT or srcEPSG must be specified")
    
        dstSR = SpatialReference()
        if dstWKT is not None:
            dstSR.ImportFromWkt(dstWKT)
        elif dstEPSG is not None:
            dstSR.ImportFromEPSG(dstEPSG)
        else:
            raise Exception("Either dstWKT or dstEPSG must be specified")
        
        transPoints = []
        
        trans = CoordinateTransformation(srcSR, dstSR)
        for pt in points:
            x, y, _ = trans.TransformPoint(pt[0], pt[1])
            transPoints.append((x, y))
        trans = None
        return transPoints
    


# ============================================================================ 

def _addUserFieldDef(self, newDataset):
        spRef = osr.SpatialReference()
        spRef.ImportFromEPSG(DEFAULT_EPSG)
        maxStrlen = LMFormat.getStrlenForDefaultOGR()
     
        newLyr = newDataset.CreateLayer('points', geom_type=ogr.wkbPoint, srs=spRef)
        if newLyr is None:
            raise LmException(JobStatus.IO_OCCURRENCE_SET_WRITE_ERROR, 
                                    'Layer creation failed')
        
        for idx, vals in self.op.columnMeta.iteritems():
            if vals is not None:
                fldname = str(vals[OccDataParser.FIELD_NAME_KEY])
                fldtype = vals[OccDataParser.FIELD_TYPE_KEY] 
                fldDef = ogr.FieldDefn(fldname, fldtype)
                if fldtype == ogr.OFTString:
                    fldDef.SetWidth(maxStrlen)
                returnVal = newLyr.CreateField(fldDef)
                if returnVal != 0:
                    raise LmException(JobStatus.IO_OCCURRENCE_SET_WRITE_ERROR, 
                                         'Failed to create field {}'.format(fldname))
                
        # Add wkt field
        fldDef = ogr.FieldDefn(LM_WKT_FIELD, ogr.OFTString)
        fldDef.SetWidth(maxStrlen)
        returnVal = newLyr.CreateField(fldDef)
        if returnVal != 0:
            raise LmException(JobStatus.IO_OCCURRENCE_SET_WRITE_ERROR, 
                                    'Failed to create field {}'.format(fldname))
        
        return newLyr
            

    # ............................................... 

def get_ccrs(filename):
    """
    Loads WKT projection string from file and return
    cartopy coordinate reference system.
    """
    inproj = osr.SpatialReference()
    proj = open(filename, 'r').readline()
    inproj.ImportFromWkt(proj)
    projcs = inproj.GetAuthorityCode('PROJCS')
    return ccrs.epsg(projcs) 

def getProjFromDict(projdict, name='', GeoCS='WGS84', convention='Proj4'):
  ''' Initialize a projected OSR SpatialReference instance from a dictionary using Proj4 conventions. 
      Valid parameters are documented here: http://trac.osgeo.org/proj/wiki/GenParms
      Projections are described here: http://www.remotesensing.org/geotiff/proj_list/ '''
  # initialize
  projection = osr.SpatialReference()
  # interpret dictionary according to convention
  if convention == 'Proj4':
    # start with projection, which is usually a string
    proj = projdict['proj']
    if proj == 'longlat' or proj == 'latlong':
      projstr = '+proj=latlong'; lproj = False
    else: 
      projstr = '+proj={0:s}'.format(proj); lproj = True 
    # loop over entries
    for key, value in projdict.items():
      if key is not 'proj':
        if not isinstance(key, str): raise TypeError
        if not isinstance(value, (float,np.inexact,int,np.integer)): raise TypeError
        # translate dict entries to string
        projstr = '{0:s} +{1:s}={2:f}'.format(projstr, key, value)
    # load projection from proj4 string
    projection.ImportFromProj4(projstr)
  else:
    raise NotImplementedError
  # more meta data
  if lproj: projection.SetProjCS(name)  # establish that this is a projected system
  # N.B.: note that the seemingly analogous method SetGeogCS() has a different function (and is not necessary)
  projection.SetWellKnownGeogCS(GeoCS)  # default reference datum/geoid
  # return finished projection object (geotransform can be inferred from coordinate vectors)
  return projection 

def get_pixel_size(self, index_tuple):
        """
        Returns X & Y sizes in metres of specified pixel as a tuple.
        N.B: Pixel ordinates are zero-based from top left
        """
        x, y = index_tuple
        logger.debug('(x, y) = (%f, %f)', x, y)
        
        native_spatial_reference = osr.SpatialReference()
        native_spatial_reference.ImportFromWkt(self.crs)
        
        latlong_spatial_reference = native_spatial_reference.CloneGeogCS()
        coord_transform_to_latlong = osr.CoordinateTransformation(native_spatial_reference, latlong_spatial_reference)

        # Determine pixel centre and edges in georeferenced coordinates
        xw = self.GeoTransform[0] + x * self.GeoTransform[1]
        yn = self.GeoTransform[3] + y * self.GeoTransform[5] 
        xc = self.GeoTransform[0] + (x + 0.5) * self.GeoTransform[1]
        yc = self.GeoTransform[3] + (y + 0.5) * self.GeoTransform[5] 
        xe = self.GeoTransform[0] + (x + 1.0) * self.GeoTransform[1]
        ys = self.GeoTransform[3] + (y + 1.0) * self.GeoTransform[5] 
        
        logger.debug('xw = %f, yn = %f, xc = %f, yc = %f, xe = %f, ys = %f', xw, yn, xc, yc, xe, ys)
        
        # Convert georeferenced coordinates to lat/lon for Vincenty
        lon1, lat1, _z = coord_transform_to_latlong.TransformPoint(xw, yc, 0)
        lon2, lat2, _z = coord_transform_to_latlong.TransformPoint(xe, yc, 0)
        logger.debug('For X size: (lon1, lat1) = (%f, %f), (lon2, lat2) = (%f, %f)', lon1, lat1, lon2, lat2)
        x_size, _az_to, _az_from = vinc_dist(earth.F, earth.A, 
                                             lat1 * RADIANS_PER_DEGREE, lon1 * RADIANS_PER_DEGREE, 
                                             lat2 * RADIANS_PER_DEGREE, lon2 * RADIANS_PER_DEGREE)
        
        lon1, lat1, _z = coord_transform_to_latlong.TransformPoint(xc, yn, 0)
        lon2, lat2, _z = coord_transform_to_latlong.TransformPoint(xc, ys, 0)
        logger.debug('For Y size: (lon1, lat1) = (%f, %f), (lon2, lat2) = (%f, %f)', lon1, lat1, lon2, lat2)
        y_size, _az_to, _az_from = vinc_dist(earth.F, earth.A, 
                                             lat1 * RADIANS_PER_DEGREE, lon1 * RADIANS_PER_DEGREE, 
                                             lat2 * RADIANS_PER_DEGREE, lon2 * RADIANS_PER_DEGREE)
        
        logger.debug('(x_size, y_size) = (%f, %f)', x_size, y_size)
        return (x_size, y_size) 

def get_spatial_ref_from_wkt(wkt_or_crs_name):
    '''
    Function to return SpatialReference object for supplied WKT
    @param wkt: Well-known text or CRS name for SpatialReference, including "EPSG:XXXX"
    @return spatial_ref: SpatialReference from WKT
    '''
    spatial_ref = SpatialReference()
    
    # Try to resolve WKT
    result = spatial_ref.ImportFromWkt(wkt_or_crs_name)
    if not result:
        return spatial_ref

    # Try to resolve CRS name - either mapped or original
    result = spatial_ref.SetWellKnownGeogCS(CRS_NAME_MAPPING.get(wkt_or_crs_name) or wkt_or_crs_name) 
    if not result:
        return spatial_ref

    # Try common formulations for UTM zones
    #TODO: Fix this so it works in the Northern hemisphere 
    modified_crs_name = re.sub('\s+', '', wkt_or_crs_name.strip().upper())
    utm_match = (re.match('(\w+)/MGAZONE(\d+)', modified_crs_name) or
                 re.match('(\w+)/(\d+)S', modified_crs_name) or
                 re.match('(EPSG:283)(\d{2})', modified_crs_name) 
                 )
    if utm_match:
        modified_crs_name = utm_match.group(1)
        utm_zone = int(utm_match.group(2))
        result = spatial_ref.SetWellKnownGeogCS(CRS_NAME_MAPPING.get(modified_crs_name) or modified_crs_name)
    if not result:
        spatial_ref.SetUTM(utm_zone, False) # Put this here to avoid potential side effects in downstream code
        return spatial_ref

    assert not result, 'Invalid WKT or CRS name' 

def get_utm_wkt(coordinate, from_wkt):
    '''
    Function to return CRS for UTM zone of specified coordinates.
    Used to transform coords to metres
    @param coordinate: single coordinate pair
    '''
    def utm_getZone(longitude):
        return (int(1 + (longitude + 180.0) / 6.0))

    def utm_isNorthern(latitude):
        if (latitude < 0.0):
            return 0
        else:
            return 1
        
    coordinate_array = np.array(coordinate).reshape((1,2))

    latlon_coord_trans = get_coordinate_transformation(
        from_wkt, 'EPSG:4326')
    latlon_coord = coordinate if latlon_coord_trans is None else latlon_coord_trans.TransformPoints(
        coordinate_array)[0][0:2]
        
    # Set UTM coordinate reference system
    utm_spatial_ref = SpatialReference()
    utm_spatial_ref.SetWellKnownGeogCS('WGS84')
    utm_spatial_ref.SetUTM(utm_getZone(
        latlon_coord[0]), utm_isNorthern(latlon_coord[1]))

    return utm_spatial_ref.ExportToWkt() 

def create_gpkg(
    gpkg_name, proj_string, size=(1, 1), geotransform=[0, 1, 0, 0, 0, -1],
    creation_options=None
):
    if os.path.exists("%s.gpkg" % gpkg_name):
        sys.stderr.write(
            "ERROR: SQLite GeoPackage '%s.gpkg' already exists.\n" % gpkg_name
        )
        sys.exit(1)

    gdal.AllRegister()
    drv = gdal.GetDriverByName("GPKG")
    try:
        gpkg = drv.Create(
            "%s.gpkg" % gpkg_name, size[0], size[1], 1, gdal.GDT_Byte,
            creation_options
        )

        proj = osr.SpatialReference()
        res = proj.SetWellKnownGeogCS(proj_string)
        if res != 0:
            if proj_string[0:4] == 'EPSG':
                proj.ImportFromEPSG(int(proj_string[5:]))
        gpkg.SetProjection(proj.ExportToWkt())
        gpkg.SetGeoTransform(geotransform)
        gpkg = None
    except Exception as e:
        sys.stderr.write(
            "ERROR: Cannot create SQLite GeoPackage '%s.gpkg'. "
            "Error message was: '%s'.\n" % (gpkg_name, e.message)
        )
        sys.exit(1) 

def proj4_virtual(self):
        if self.wkt_virtual is None:
            return None # pragma: no cover
        return osr.SpatialReference(self.wkt_virtual).ExportToProj4() 

def proj4_stored(self):
        if self.wkt_stored is None:
            return None # pragma: no cover
        return osr.SpatialReference(self.wkt_stored).ExportToProj4() 

def wkt_same(a, b):
    """Are two wkt equivalent"""
    if a == b:
        return True
    sra = osr.SpatialReference(a)
    srb = osr.SpatialReference(b)
    return bool(sra.IsSame(srb)) 

def sreq(*items, verbose_on=None):
    """SR items are all equal"""
    v = tuple(map(int, gdal.__version__.split('.')))
    for a, b in itertools.combinations(items, 2):
        a = osr.SpatialReference(osr.GetUserInputAsWKT(a))
        if v < (3,):
            a.StripCTParms()
        a = a.ExportToProj4()
        a = osr.SpatialReference(osr.GetUserInputAsWKT(a))
        if v < (3,):
            a.StripCTParms()

        b = osr.SpatialReference(osr.GetUserInputAsWKT(b))
        if v < (3,):
            b.StripCTParms()
        b = b.ExportToProj4()
        b = osr.SpatialReference(osr.GetUserInputAsWKT(b))
        if v < (3,):
            b.StripCTParms()

        res = bool(a.IsSame(b))
        if res is verbose_on:
            print('')
            print(a.ExportToPrettyWkt())
            print('---------- vs ----------')
            print(b.ExportToPrettyWkt())
            print('')

        if not res:
            return False
    return True 

def proj4(self):
        """Dataset's work spatial reference in WKT proj4.
        Returns None if `mode 1`.
        """
        if self._back.wkt_work is None:
            return None
        return osr.SpatialReference(self._back.wkt_work).ExportToProj4() 

def __init__(self, file, band=1):
        self.dataset = gdal.Open(file)
        self.geotransform = self.dataset.GetGeoTransform()
        
        srs = osr.SpatialReference()
        srs.ImportFromWkt(self.dataset.GetProjection())
        self.proj = pyproj.Proj(srs.ExportToProj4())
        
        band = self.dataset.GetRasterBand(band) 

def GetGeoInfo(FileName):
    """This gets information from the raster file using gdal

    Args:
        FileName (str): The filename (with path and extension) of the raster.

    Return:
        float: A vector that contains:
            * NDV: the nodata values
            * xsize: cellsize in x direction
            * ysize: cellsize in y direction
            * GeoT: the tranform (a string)
            * Projection: the Projection (a string)
            * DataType: The type of data (an int explaing the bits of each data element)

    Author: SMM
    """


    if exists(FileName) is False:
            raise Exception('[Errno 2] No such file or directory: \'' + FileName + '\'')


    SourceDS = gdal.Open(FileName, gdal.GA_ReadOnly)
    if SourceDS == None:
        raise Exception("Unable to read the data file")

    NDV = SourceDS.GetRasterBand(1).GetNoDataValue()
    xsize = SourceDS.RasterXSize
    ysize = SourceDS.RasterYSize
    GeoT = SourceDS.GetGeoTransform()
    Projection = osr.SpatialReference()
    Projection.ImportFromWkt(SourceDS.GetProjectionRef())
    DataType = SourceDS.GetRasterBand(1).DataType
    DataType = gdal.GetDataTypeName(DataType)

    return NDV, xsize, ysize, GeoT, Projection, DataType
#==============================================================================

#==============================================================================
# This gets the UTM zone, if it exists 

def array2raster(rasterfn,newRasterfn,array,driver_name = "ENVI", noDataValue = -9999):
    """Takes an array and writes to a GDAL compatible raster. It needs another raster to map the dimensions.

    Args:
        FileName (str): The filename (with path and extension) of a raster that has the same dimensions as the raster to be written.
        newRasterfn (str): The filename (with path and extension) of the new raster.
        array (np.array): The array to be written
        driver_name (str): The type of raster to write. Default is ENVI since that is the LSDTOpoTools format
        noDataValue (float): The no data value

    Return:
        np.array: A numpy array with the data from the raster.

    Author: SMM
    """

    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    cols = raster.RasterXSize
    rows = raster.RasterYSize

    driver = gdal.GetDriverByName(driver_name)
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outRaster.GetRasterBand(1).SetNoDataValue( noDataValue )
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()
#============================================================================== 

def read(self, infile):
		# never want to change the original input raster so use read only constant
		self.infile = infile
		self.raster = gdal.Open(self.infile, GA_ReadOnly)
		self.band = self.raster.GetRasterBand(1)
		self.NDV = self.band.GetNoDataValue()
		self.x = self.band.XSize
		self.y = self.band.YSize
		DataType = self.band.DataType
		self.DataType = gdal.GetDataTypeName(DataType)
		self.GeoT = self.raster.GetGeoTransform()
		# Problem
		self.prj = osr.SpatialReference()
		self.prj.ImportFromWkt(self.raster.GetProjectionRef()) 

def polygonize(img,shp_path):
    # mapping between gdal type and ogr field type
    type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                    gdal.GDT_UInt16: ogr.OFTInteger,
                    gdal.GDT_Int16: ogr.OFTInteger,
                    gdal.GDT_UInt32: ogr.OFTInteger,
                    gdal.GDT_Int32: ogr.OFTInteger,
                    gdal.GDT_Float32: ogr.OFTReal,
                    gdal.GDT_Float64: ogr.OFTReal,
                    gdal.GDT_CInt16: ogr.OFTInteger,
                    gdal.GDT_CInt32: ogr.OFTInteger,
                    gdal.GDT_CFloat32: ogr.OFTReal,
                    gdal.GDT_CFloat64: ogr.OFTReal}

    ds = gdal.Open(img)
    prj = ds.GetProjection()
    srcband = ds.GetRasterBand(1)

    dst_layername = "Shape"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference(wkt=prj)

    dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
    # raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType])
    raster_field = ogr.FieldDefn('id', type_mapping[gdal.GDT_Int32])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# extract sinks from dem 

def wgs84_to_modissinu ( lon, lat ):
	wgs84 = osr.SpatialReference( )  # Define a SpatialReference object
	wgs84.ImportFromEPSG( 4326 )
	# Define the modis data projection, sinusoidal #sinu_modis.ImportFromEPSG( 27700)
	modis_sinu = osr.SpatialReference() 
	modis_sinu.ImportFromProj4 ( "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 " + \
                      "+a=6371007.181 +b=6371007.181 " + \
                      "+units=m +no_defs" )
	tx = osr.CoordinateTransformation(wgs84, modis_sinu)
	modis_x, modis_y, modis_z = tx.TransformPoint ( lon, lat )
	return modis_x, modis_y, modis_z 

def modissinu_to_wgs84 ( modis_x, modis_y ):
	from osgeo import osr
	wgs84 = osr.SpatialReference( )
	wgs84.ImportFromEPSG( 4326 )
	modis_sinu = osr.SpatialReference() 
	modis_sinu.ImportFromProj4 ( "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 " + \
                      "+a=6371007.181 +b=6371007.181 " + \
                      "+units=m +no_defs" )
	tx = osr.CoordinateTransformation(modis_sinu, wgs84)
	lon, lat, modis_z = tx.TransformPoint ( modis_x, modis_y ) 

def wgs84_to_modissinu ( lon, lat ):
	wgs84 = osr.SpatialReference( )
	wgs84.ImportFromEPSG( 4326 )
	modis_sinu = osr.SpatialReference() 
	modis_sinu.ImportFromProj4 ( "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 " + \
                      "+a=6371007.181 +b=6371007.181 " + \
                      "+units=m +no_defs" )
	tx = osr.CoordinateTransformation(wgs84, modis_sinu)
	modis_x, modis_y, modis_z = tx.TransformPoint ( lon, lat )
	return modis_x, modis_y, modis_z 

def modissinu_to_wgs84 ( modis_x, modis_y ):
	from osgeo import osr
	wgs84 = osr.SpatialReference( )
	wgs84.ImportFromEPSG( 4326 )
	modis_sinu = osr.SpatialReference() 
	modis_sinu.ImportFromProj4 ( "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 " + \
                      "+a=6371007.181 +b=6371007.181 " + \
                      "+units=m +no_defs" )
	tx = osr.CoordinateTransformation(modis_sinu, wgs84)
	lon, lat, modis_z = tx.TransformPoint ( modis_x, modis_y ) 

def output_file(output_name,output_array):
  
    driver = gdal.GetDriverByName("GTiff")
    outDataset = driver.Create(output_name, ncols, nrows, 1, gdal.GDT_Float32)
    outDataset.SetGeoTransform(geo_transform_1)
    #srs = osr.SpatialReference()
    #srs.ImportFromEPSG(4326) #WGS84 lat long.
    outDataset.SetProjection(proj_1)
    
    outBand = outDataset.GetRasterBand(1)
    outBand.WriteArray(output_array,0,0)
    
    outBand.SetNoDataValue(np.nan)
    outBand.FlushCache() 

def preprocess_sen2_images(l2_dir, out_dir, l1_dir, cloud_threshold=60, buffer_size=0, epsg=None,
                           bands=("B08", "B04", "B03", "B02"), out_resolution=10):
    """For every .SAFE folder in in_dir, stacks band 2,3,4 and 8  bands into a single geotif, creates a cloudmask from
    the combined fmask and sen2cor cloudmasks and reprojects to a given EPSG if provided"""
    safe_file_path_list = [os.path.join(l2_dir, safe_file_path) for safe_file_path in os.listdir(l2_dir)]
    for l2_safe_file in safe_file_path_list:
        with TemporaryDirectory() as temp_dir:
            log.info("----------------------------------------------------")
            log.info("Merging 10m bands in SAFE dir: {}".format(l2_safe_file))
            temp_path = os.path.join(temp_dir, get_sen_2_granule_id(l2_safe_file)) + ".tif"
            log.info("Output file: {}".format(temp_path))
            stack_sentinel_2_bands(l2_safe_file, temp_path, bands=bands, out_resolution=out_resolution)

            log.info("Creating cloudmask for {}".format(temp_path))
            l1_safe_file = get_l1_safe_file(l2_safe_file, l1_dir)
            mask_path = get_mask_path(temp_path)
            create_mask_from_sen2cor_and_fmask(l1_safe_file, l2_safe_file, mask_path, buffer_size=buffer_size)
            log.info("Cloudmask created")

            out_path = os.path.join(out_dir, os.path.basename(temp_path))
            out_mask_path = os.path.join(out_dir, os.path.basename(mask_path))

            if epsg:
                log.info("Reprojecting images to {}".format(epsg))
                proj = osr.SpatialReference()
                proj.ImportFromEPSG(epsg)
                wkt = proj.ExportToWkt()
                reproject_image(temp_path, out_path, wkt)
                reproject_image(mask_path, out_mask_path, wkt)
            else:
                log.info("Moving images to {}".format(out_dir))
                shutil.move(temp_path, out_path)
                shutil.move(mask_path, out_mask_path) 

def reproject_geotransform(in_gt, old_proj_wkt, new_proj_wkt):
    """
    Reprojects a geotransform from the old projection to a new projection. See
    [https://gdal.org/user/raster_data_model.html]

    Parameters
    ----------
    in_gt
        A six-element numpy array, usually an output from gdal_image.GetGeoTransform()
    old_proj_wkt
        The projection of the old geotransform in well-known text.
    new_proj_wkt
        The projection of the new geotrasform in well-known text.
    Returns
    -------
    out_gt
        The geotransform in the new projection

    """
    old_proj = osr.SpatialReference()
    new_proj = osr.SpatialReference()
    old_proj.ImportFromWkt(old_proj_wkt)
    new_proj.ImportFromWkt(new_proj_wkt)
    transform = osr.CoordinateTransformation(old_proj, new_proj)
    (ulx, uly, _) = transform.TransformPoint(in_gt[0], in_gt[3])
    out_gt = (ulx, in_gt[1], in_gt[2], uly, in_gt[4], in_gt[5])
    return out_gt 

def write_geometry(geometry, out_path, srs_id=4326):
    """
    Saves the geometry in an ogr.Geometry object to a shapefile.

    Parameters
    ----------
    geometry
        An ogr.Geometry object
    out_path
        The location to save the output shapefile
    srs_id
        The projection of the output shapefile. Can be an EPSG number or a WKT string.

    Notes
    -----
    The shapefile consists of one layer named 'geometry'.


    """
    # TODO: Fix this needing an extra filepath on the end
    driver = ogr.GetDriverByName("ESRI Shapefile")
    data_source = driver.CreateDataSource(out_path)
    srs = osr.SpatialReference()
    if type(srs_id) is int:
        srs.ImportFromEPSG(srs_id)
    if type(srs_id) is str:
        srs.ImportFromWkt(srs_id)
    layer = data_source.CreateLayer(
        "geometry",
        srs,
        geom_type=geometry.GetGeometryType())
    feature_def = layer.GetLayerDefn()
    feature = ogr.Feature(feature_def)
    feature.SetGeometry(geometry)
    layer.CreateFeature(feature)
    data_source.FlushCache()
    data_source = None 

def __str__(self):
        """
        returns the projection spec.
        """
        pszProjection = self.data_set.GetProjectionRef()
        if pszProjection is not None:
            hSRS = osr.SpatialReference()
            if hSRS.ImportFromWkt(pszProjection ) == gdal.CE_None:
                pszPrettyWkt = hSRS.ExportToPrettyWkt(False)
                return pszPrettyWkt
            else:
                return pszProjection 

def srs(self):
        """SpatialReference: Spatial reference system used by the loaded raster file."""
        if self._srs is None:
            srs = osr.SpatialReference()
            srs.ImportFromWkt(self._ds.GetProjection())
            self._srs = srs
        return self._srs 

def srs_handler(ctx, param, value):
    out_srs = osr.SpatialReference()
    out_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    if out_srs.SetFromUserInput(value) != 0:
        raise ValueError("Failed to process SRS definition: %s" % value)
    return out_srs 

def test_maskedrasterreader(classraster_filepath):
    reader = MaskedRasterReader(classraster_filepath)
    assert reader
    assert isinstance(reader.srs, osr.SpatialReference)
    assert reader.bbox == Bbox(727000.0, 6171000.0, 728000.0, 6172000.0)
    pix_win = reader.bbox_to_pixel_window(reader.bbox)
    assert pix_win == (0, 0, 500, 500)
    pix_win = reader.bbox_to_pixel_window(
        Bbox(727500.0, 6171600.0, 727600.0, 6171750.0)
    )
    assert pix_win == (250, 125, 50, 75)
    # Try using a point
    pnt = ogr.Geometry(ogr.wkbPoint)
    pnt.AddPoint(727500.0, 6171600.0)
    pnt.AssignSpatialReference(reader.srs)
    data = reader.read_2d(pnt)
    np.testing.assert_array_equal(data, np.ma.empty(shape=(0, 0)))
    # Now use a proper polygon
    poly = pnt.Buffer(100)
    data = reader.read_2d(poly)
    assert data.shape == (100, 100)
    # Check that we have a mask
    assert int(np.sum(data.mask)) == 2140
    # Check data without mask
    assert int(np.sum(data.data)) == 25412
    # Check data with mask. (Must be less than unmasked)
    assert int(np.sum(data.compressed())) == 20432

    # Test read_flattened
    flat_data = reader.read_flattened(poly)
    # read_2d returns a 100x100 of which 2140 are outside the poly
    assert flat_data.shape == (100 * 100 - 2140,)
    assert np.ma.is_masked(flat_data)
    assert int(np.ma.sum(flat_data)) == 20432 

def __init__(self, config):
    default_config = {
      "projection": "merc",
      "longitude0": 0
    }
    default_config.update(config)
    self.config = default_config

    self.spatialRef = osr.SpatialReference()
    projString = '+proj='+str(self.config['projection'])+' +a=6381372 +b=6381372 +lat_0=0'
    #if 'emulate_longitude0' in self.config and not self.config['emulate_longitude0']:
    projString += ' +lon_0='+str(self.config['longitude0'])
    self.spatialRef.ImportFromProj4(projString) 

def update_map(mapname, newfilepath, get_proj=None):
	if mapname in maps_paths and maps_paths[mapname] == newfilepath:
		return
	dataset = gdal.Open(newfilepath)
	if dataset:
		maps[mapname] = dataset
		maps_paths[mapname] = newfilepath
		proj = osr.SpatialReference()
		proj.ImportFromWkt(dataset.GetProjection())
		if proj.Validate() != 0 and get_proj:
			epsg = get_proj(mapname)
			proj.ImportFromEPSG(epsg)
			dataset.SetProjection(proj.ExportToWkt())