Python osgeo.osr.SpatialReference() Examples

def convert_pixgwktList_to_wgs84wktList(inputRaster, wktPolygonPixList):
    ## returns # [[GeoWKT, PixWKT], ...]
    wgs84WKTList=[]
    if os.path.isfile(inputRaster):
        srcRaster = gdal.Open(inputRaster)
        targetSR = osr.SpatialReference()
        targetSR.ImportFromWkt(srcRaster.GetProjectionRef())
        geomTransform = srcRaster.GetGeoTransform()

    for wktPolygonPix in wktPolygonPixList:
        geom_wkt_list = pixelWKTToGeoWKT(wktPolygonPix, inputRaster, targetSR='',
                                         geomTransform=geomTransform,
                                         breakMultiPolygonPix=False)

        wgs84WKTList.extend(geom_wkt_list)

    # [[GeoWKT, PixWKT], ...]
    return wgs84WKTList 

def polygonize(self,rasterTemp,outShp):
            
            sourceRaster = gdal.Open(rasterTemp)
            band = sourceRaster.GetRasterBand(1)
            driver = ogr.GetDriverByName("ESRI Shapefile")
            # If shapefile already exist, delete it
            if os.path.exists(outShp):
                driver.DeleteDataSource(outShp)
                
            outDatasource = driver.CreateDataSource(outShp)            
            # get proj from raster            
            srs = osr.SpatialReference()
            srs.ImportFromWkt( sourceRaster.GetProjectionRef() )
            # create layer with proj
            outLayer = outDatasource.CreateLayer(outShp,srs)
            # Add class column (1,2...) to shapefile
      
            newField = ogr.FieldDefn('Class', ogr.OFTInteger)
            outLayer.CreateField(newField)
            
            gdal.Polygonize(band, None,outLayer, 0,[],callback=None)  
            
            outDatasource.Destroy()
            sourceRaster=None
            band=None
                    
            
            ioShpFile = ogr.Open(outShp, update = 1)
            
            
            lyr = ioShpFile.GetLayerByIndex(0)
            
            lyr.ResetReading()    

            for i in lyr:
                lyr.SetFeature(i)
            # if area is less than inMinSize or if it isn't forest, remove polygon 
                if i.GetField('Class')!=1:
                    lyr.DeleteFeature(i.GetFID())        
            ioShpFile.Destroy()
                        
            return outShp 

def toWgs(mgrs):
    """ Converts an MGRS coordinate string to geodetic (latitude and longitude)
    coordinates

    @param mgrs - MGRS coordinate string
    @returns - tuple containning latitude and longitude values
    """
    if _checkZone(mgrs):
        zone, hemisphere, easting, northing = _mgrsToUtm(mgrs)
    else:
        zone, hemisphere, easting, northing = _mgrsToUps(mgrs)

    epsg = _epsgForUtm(zone, hemisphere)
    src = osr.SpatialReference()
    src.ImportFromEPSG(epsg)
    dst = osr.SpatialReference()
    dst.ImportFromEPSG(4326)
    ct = osr.CoordinateTransformation(src, dst)
    longitude, latitude, z = ct.TransformPoint(easting, northing)

    return latitude, longitude 

def raster2array(rasterfn):
    raster = gdal.Open(rasterfn)
    band = raster.GetRasterBand(1)
    nodata = band.GetNoDataValue()
    array = band.ReadAsArray()
    
    proj = raster.GetProjection()
    inproj = osr.SpatialReference()
    inproj.ImportFromWkt(proj)
    
    geoTransform = raster.GetGeoTransform()
    minx = geoTransform[0]
    maxy = geoTransform[3]
    maxx = minx + geoTransform[1]*raster.RasterXSize
    miny = maxy + geoTransform[5]*raster.RasterYSize
    extent =  [minx, maxx, miny, maxy]
    pixelSizeXY = [geoTransform[1], geoTransform[5]]
    del raster, band
    return [array, nodata, extent, inproj, pixelSizeXY]

#clip a raster by vector 

def _get_raster_attributes(self,ds,attrs):
        """Collect raster attributes
        """

        attrs["type"] = "raster"

        geotransform = ds.GetGeoTransform()
        attrs["extent"] = (geotransform[0], 
                           geotransform[3]+(geotransform[5]*ds.RasterYSize), 
                           geotransform[0]+(geotransform[1]*ds.RasterXSize),
                           geotransform[3])

        sr = osr.SpatialReference()
        sr.ImportFromWkt(ds.GetProjectionRef())
        attrs["prj"] = self._get_prj(sr)

        attrs["features_count"] = "%s raster, %dx%d cells" % \
                (ds.RasterCount, ds.RasterXSize, ds.RasterYSize)

        return attrs 

def read_gdal_projection(dset):
    """Get a projection (OSR object) from a GDAL dataset.

    Parameters
    ----------
    dset : gdal.Dataset

    Returns
    -------
    srs : OSR.SpatialReference
        dataset projection object

    Examples
    --------

    See :ref:`notebooks/classify/wradlib_clutter_cloud_example.ipynb`.

    """
    wkt = dset.GetProjection()
    srs = osr.SpatialReference()
    srs.ImportFromWkt(wkt)
    # src = None
    return srs 

def proj4_to_osr(proj4str):
    """Transform a proj4 string to an osr spatial reference object

    Parameters
    ----------
    proj4str : string
        Proj4 string describing projection

    Examples
    --------

    See :ref:`notebooks/radolan/radolan_grid.ipynb#PROJ.4`.

    """
    proj = None
    if proj4str:
        proj = osr.SpatialReference()
        proj.ImportFromProj4(proj4str)
    else:
        proj = get_default_projection()
    return proj 

def epsg_to_osr(epsg=None):
    """Create osr spatial reference object from EPSG number

    .. versionadded:: 0.6.0

    Parameters
    ----------
    epsg : int
        EPSG-Number defining the coordinate system

    Returns
    -------
    proj : osr.SpatialReference
        GDAL/OSR object defining projection
    """
    proj = None
    if epsg:
        proj = osr.SpatialReference()
        proj.ImportFromEPSG(epsg)
    else:
        proj = get_default_projection()
    return proj 

def wkt_to_osr(wkt=None):
    """Create osr spatial reference object from WKT string

    .. versionadded:: 0.6.0

    Parameters
    ----------
    wkt : str
        WTK string defining the coordinate reference system

    Returns
    -------
    proj : osr.SpatialReference
        GDAL/OSR object defining projection

    """
    proj = None
    if wkt:
        proj = osr.SpatialReference()
        proj.ImportFromWkt(wkt)
    else:
        proj = get_default_projection()
    return proj 

def test_get_radolan_grid_equality(self):
        # create radolan projection osr object
        scale = (1. + np.sin(np.radians(60.))) / (1. + np.sin(np.radians(90.)))
        dwd_string = ('+proj=stere +lat_0=90 +lat_ts=90 +lon_0=10 '
                      '+k={0:10.8f} +x_0=0 +y_0=0 +a=6370040 +b=6370040 '
                      '+to_meter=1000 +no_defs'.format(scale))
        proj_stereo = georef.proj4_to_osr(dwd_string)

        # create wgs84 projection osr object
        proj_wgs = osr.SpatialReference()
        proj_wgs.ImportFromEPSG(4326)

        # transform radolan polar stereographic projection to wgs84 and wgs84
        # to polar stereographic
        # using osr transformation routines
        radolan_grid_ll = georef.reproject(self.radolan_grid_xy,
                                           projection_source=proj_stereo,
                                           projection_target=proj_wgs)
        radolan_grid_xy = georef.reproject(self.radolan_grid_ll,
                                           projection_source=proj_wgs,
                                           projection_target=proj_stereo)

        # check source and target arrays for equality
        self.assertTrue(np.allclose(radolan_grid_ll, self.radolan_grid_ll))
        self.assertTrue(np.allclose(radolan_grid_xy, self.radolan_grid_xy)) 

def addMetricDistanceToEdge(x1,y1,x2,y2, epsgCode):
	# we assume initial epsg is wsg84 (merctor projection)
	
	if epsgCode != 4326:
		sourceProjection = osr.SpatialReference()
		sourceProjection.ImportFromEPSG(4326)
		destinationProjection = osr.SpatialReference()
		destinationProjection.ImportFromEPSG(epsgCode) # https://epsg.io/2154
		coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	
	line = ogr.Geometry(ogr.wkbLineString)
	line.AddPoint(x1,y1)
	line.AddPoint(x2,y2)
	if epsgCode != 4326:
		line.Transform(coordTrans)
	length = line.Length()

	return length 

def addMetricDistanceToEdges(graph, epsgCode):
	# we assume initial epsg is wsg84 (merctor projection)
	metricDistance = {}
	
	if epsgCode != 4326:
		sourceProjection = osr.SpatialReference()
		sourceProjection.ImportFromEPSG(4326)
		destinationProjection = osr.SpatialReference()
		destinationProjection.ImportFromEPSG(epsgCode) # https://epsg.io/2154
		coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	
	for edge in graph.edges():
		node1, node2 =  edge
		line = ogr.Geometry(ogr.wkbLineString)
		line.AddPoint(graph.node[node1]['longitude'], graph.node[node1]['latitude'])
		line.AddPoint(graph.node[node2]['longitude'], graph.node[node2]['latitude'])
		if epsgCode != 4326:
			line.Transform(coordTrans)
		length = line.Length()
		metricDistance[edge] = length
	nx.set_edge_attributes(graph, 'length', metricDistance)

	return graph 

def addMetricDistanceToEdges(graph, epsgCode):
	# we assume initial epsg is wsg84 (merctor projection)
	metricDistance = {}
	
	if epsgCode != 4326:
		sourceProjection = osr.SpatialReference()
		sourceProjection.ImportFromEPSG(4326)
		destinationProjection = osr.SpatialReference()
		destinationProjection.ImportFromEPSG(epsgCode) # https://epsg.io/2154
		coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	
	for edge in graph.edges():
		node1, node2 =  edge
		line = ogr.Geometry(ogr.wkbLineString)
		line.AddPoint(graph.node[node1]['longitude'], graph.node[node1]['latitude'])
		line.AddPoint(graph.node[node2]['longitude'], graph.node[node2]['latitude'])
		if epsgCode != 4326:
			line.Transform(coordTrans)
		length = line.Length()
		metricDistance[edge] = length
	nx.set_edge_attributes(graph, 'length', metricDistance)

	return graph 

def toWgs(mgrs):
    """ Converts an MGRS coordinate string to geodetic (latitude and longitude)
    coordinates

    @param mgrs - MGRS coordinate string
    @returns - tuple containning latitude and longitude values
    """
    if _checkZone(mgrs):
        zone, hemisphere, easting, northing = _mgrsToUtm(mgrs)
    else:
        zone, hemisphere, easting, northing = _mgrsToUps(mgrs)

    epsg = _epsgForUtm(zone, hemisphere)
    src = osr.SpatialReference()
    src.ImportFromEPSG(epsg)
    dst = osr.SpatialReference()
    dst.ImportFromEPSG(4326)
    ct = osr.CoordinateTransformation(src, dst)
    longitude, latitude, z = ct.TransformPoint(easting, northing)

    return latitude, longitude 

def array2raster(newRasterFilename,rasterOrigin,pixelWidth,pixelHeight,array,dataType):
	cols=array.shape[1]
	rows=array.shape[0]
	originX=rasterOrigin[0]
	originY=rasterOrigin[1]
	driver=gdal.GetDriverByName('GTiff')
	outRaster = driver.Create(newRasterFilename, cols, rows, 1, dataType)
	outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
	outband = outRaster.GetRasterBand(1)
	outband.WriteArray(array)
	outRasterSRS = osr.SpatialReference()
	outRasterSRS.ImportFromEPSG(32645)# this is the EPSG code for Nepal, should be changed for other locations
	outRaster.SetProjection(outRasterSRS.ExportToWkt())
	outband.FlushCache()

#takes a series of rasters and clips them to minExtent, then returns them as numpy arrays 

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 srs(self):
        sr = SpatialReference()
        sr.ImportFromEPSG(self.srsid)
        return sr.ExportToWkt() 

def get_bbox(shapefile):
    """ Gets the bounding box of a shapefile in EPSG 4326.
        If shapefile is not in WGS84, bounds are reprojected.
    """
    driver = ogr.GetDriverByName('ESRI Shapefile')
    # 0 means read, 1 means write
    data_source = driver.Open(shapefile, 0)
    # Check to see if shapefile is found.
    if data_source is None:
        failure('Could not open %s' % (shapefile))
    else:
        layer = data_source.GetLayer()
        shape_bbox = layer.GetExtent()

        spatialRef = layer.GetSpatialRef()
        target = osr.SpatialReference()
        target.ImportFromEPSG(4326)

        # this check for non-WGS84 projections gets some false positives, but that's okay
        if target.ExportToProj4() != spatialRef.ExportToProj4():
            transform = osr.CoordinateTransformation(spatialRef, target)
            point1 = ogr.Geometry(ogr.wkbPoint)
            point1.AddPoint(shape_bbox[0], shape_bbox[2])
            point2 = ogr.Geometry(ogr.wkbPoint)
            point2.AddPoint(shape_bbox[1], shape_bbox[3])
            point1.Transform(transform)
            point2.Transform(transform)
            shape_bbox = [point1.GetX(), point2.GetX(), point1.GetY(), point2.GetY()]

        return shape_bbox 

def lat_lon_to_pixel(raster_dataset, location):
    """From zacharybears.com/using-python-to-translate-latlon-locations-to-pixels-on-a-geotiff/."""
    ds = raster_dataset
    gt = ds.GetGeoTransform()
    srs = osr.SpatialReference()
    srs.ImportFromWkt(ds.GetProjection())
    srs_lat_lon = srs.CloneGeogCS()
    ct = osr.CoordinateTransformation(srs_lat_lon, srs)
    new_location = [None, None]
    # Change the point locations into the GeoTransform space
    (new_location[1], new_location[0], holder) = ct.TransformPoint(location[1], location[0])
    # Translate the x and y coordinates into pixel values
    x = (new_location[1] - gt[0]) / gt[1]
    y = (new_location[0] - gt[3]) / gt[5]
    return (int(x), int(y)) 

def pixel_to_lat_lon(raster_dataset, col, row):
    """From zacharybears.com/using-python-to-translate-latlon-locations-to-pixels-on-a-geotiff/."""
    ds = raster_dataset
    gt = ds.GetGeoTransform()
    srs = osr.SpatialReference()
    srs.ImportFromWkt(ds.GetProjection())
    srs_lat_lon = srs.CloneGeogCS()
    ct = osr.CoordinateTransformation(srs, srs_lat_lon)
    ulon = col * gt[1] + gt[0]
    ulat = row * gt[5] + gt[3]
    # Transform the point into the GeoTransform space
    (lon, lat, holder) = ct.TransformPoint(ulon, ulat)
    return (lat, lon) 

def pixel_to_lat_lon_web_mercator(raster_dataset, col, row):
    """Convert a pixel on the raster_dataset to web mercator (epsg:3857)."""
    spatial_reference = osr.SpatialReference()
    spatial_reference.ImportFromWkt(raster_dataset.GetProjection())
    ds_spatial_reference_proj_string = spatial_reference.ExportToProj4()
    in_proj = Proj(ds_spatial_reference_proj_string)
    out_proj = Proj(init='epsg:3857')

    geo_transform = raster_dataset.GetGeoTransform()
    ulon = col * geo_transform[1] + geo_transform[0]
    ulat = row * geo_transform[5] + geo_transform[3]

    x2, y2 = transform(in_proj, out_proj, ulon, ulat)
    x2, y2 = out_proj(x2, y2, inverse=True)
    return x2, y2 

def buildTindex(rasterFolder, rasterExtention='.tif'):
    rasterList = glob.glob(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))
    print(rasterList)

    print(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))

    memDriver = ogr.GetDriverByName('MEMORY')
    gTindex = memDriver.CreateDataSource('gTindex')
    srcImage = gdal.Open(rasterList[0])
    spat_ref = osr.SpatialReference()
    spat_ref.SetProjection(srcImage.GetProjection())
    gTindexLayer = gTindex.CreateLayer("gtindexlayer", spat_ref, geom_type=ogr.wkbPolygon)

    # Add an ID field
    idField = ogr.FieldDefn("location", ogr.OFTString)
    gTindexLayer.CreateField(idField)

    # Create the feature and set values
    featureDefn = gTindexLayer.GetLayerDefn()



    for rasterFile in rasterList:
        srcImage = gdal.Open(rasterFile)

        geoTrans, polyToCut, ulX, ulY, lrX, lrY = gT.getRasterExtent(srcImage)

        feature = ogr.Feature(featureDefn)
        feature.SetGeometry(polyToCut)
        feature.SetField("location", rasterFile)
        gTindexLayer.CreateFeature(feature)
        feature = None


    return gTindex, gTindexLayer 

def latlon2pixel(lat, lon, input_raster='', targetsr='', geom_transform=''):
    # type: (object, object, object, object, object) -> object

    sourcesr = osr.SpatialReference()
    sourcesr.ImportFromEPSG(4326)

    geom = ogr.Geometry(ogr.wkbPoint)
    geom.AddPoint(lon, lat)

    if targetsr == '':
        src_raster = gdal.Open(input_raster)
        targetsr = osr.SpatialReference()
        targetsr.ImportFromWkt(src_raster.GetProjectionRef())
    coord_trans = osr.CoordinateTransformation(sourcesr, targetsr)
    if geom_transform == '':
        src_raster = gdal.Open(input_raster)
        transform = src_raster.GetGeoTransform()
    else:
        transform = geom_transform

    x_origin = transform[0]
    # print(x_origin)
    y_origin = transform[3]
    # print(y_origin)
    pixel_width = transform[1]
    # print(pixel_width)
    pixel_height = transform[5]
    # print(pixel_height)
    geom.Transform(coord_trans)
    # print(geom.GetPoint())
    x_pix = (geom.GetPoint()[0] - x_origin) / pixel_width
    y_pix = (geom.GetPoint()[1] - y_origin) / pixel_height

    return (x_pix, y_pix) 

def returnBoundBox(xOff, yOff, pixDim, inputRaster, targetSR='', pixelSpace=False):
    # Returns Polygon for a specific square defined by a center Pixel and
    # number of pixels in each dimension.
    # Leave targetSR as empty string '' or specify it as a osr.SpatialReference()
    # targetSR = osr.SpatialReference()
    # targetSR.ImportFromEPSG(4326)
    if targetSR == '':
        targetSR = osr.SpatialReference()
        targetSR.ImportFromEPSG(4326)
    xCord = [xOff - pixDim / 2, xOff - pixDim / 2, xOff + pixDim / 2,
             xOff + pixDim / 2, xOff - pixDim / 2]
    yCord = [yOff - pixDim / 2, yOff + pixDim / 2, yOff + pixDim / 2,
             yOff - pixDim / 2, yOff - pixDim / 2]

    ring = ogr.Geometry(ogr.wkbLinearRing)
    for idx in xrange(len(xCord)):
        if pixelSpace == False:
            geom = pixelToGeoCoord(xCord[idx], yCord[idx], inputRaster)
            ring.AddPoint(geom[0], geom[1], 0)
        else:
            ring.AddPoint(xCord[idx], yCord[idx], 0)

    poly = ogr.Geometry(ogr.wkbPolygon)
    if pixelSpace == False:
        poly.AssignSpatialReference(targetSR)

    poly.AddGeometry(ring)

    return poly 

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 toMgrs(latitude, longitude, precision=5):
    """ Converts geodetic (latitude and longitude) coordinates to an MGRS
    coordinate string, according to the current ellipsoid parameters.

    @param latitude - latitude value
    @param longitude - longitude value
    @param precision - precision level of MGRS string
    @returns - MGRS coordinate string
    """
    if math.fabs(latitude) > 90:
        raise MgrsException('Latitude outside of valid range (-90 to 90 degrees).')

    if (longitude < -180) or (longitude > 360):
        raise MgrsException('Longitude outside of valid range (-180 to 360 degrees).')

    if (precision < 0) or (precision > MAX_PRECISION):
        raise MgrsException('The precision must be between 0 and 5 inclusive.')

    hemisphere, zone, epsg = _epsgForWgs(latitude, longitude)
    src = osr.SpatialReference()
    src.ImportFromEPSG(4326)
    dst = osr.SpatialReference()
    dst.ImportFromEPSG(epsg)
    ct = osr.CoordinateTransformation(src, dst)
    x, y, z = ct.TransformPoint(longitude, latitude)

    if (latitude < -80) or (latitude > 84):
        # Convert to UPS
        mgrs = _upsToMgrs(hemisphere, x, y, precision)
    else:
        # Convert to UTM
        mgrs = _utmToMgrs(zone, hemisphere, latitude, longitude, x, y, precision)

    return mgrs 

def latlng_bounding_box_from_ds(ds):
    """
    Given a GDAL dataset, computes its bounding box int lat/lng
    """
    wgs84_sr = osr.SpatialReference()
    wgs84_sr.ImportFromEPSG(4326)

    ds_sr = osr.SpatialReference()
    ds_sr.ImportFromWkt(ds.GetProjection())

    ds_to_wgs84 = osr.CoordinateTransformation(ds_sr, wgs84_sr)

    geot = ds.GetGeoTransform()

    def _xy2latlng(x, y):
        x_geo = geot[0] + x * geot[1] + y * geot[2]
        y_geo = geot[3] + x * geot[4] + y * geot[5]
        lng, lat, _ = ds_to_wgs84.TransformPoint(x_geo, y_geo)
        return lat, lng
    poly = [
        _xy2latlng(0, 0),
        _xy2latlng(ds.RasterXSize, 0),
        _xy2latlng(ds.RasterXSize, ds.RasterYSize),
        _xy2latlng(0, ds.RasterYSize)
    ]
    return poly 

def latlng_bbox_from_ds(ds):
    wgs84_sr = osr.SpatialReference()
    wgs84_sr.ImportFromEPSG(4326)

    ds_sr = osr.SpatialReference()
    ds_sr.ImportFromWkt(ds.GetProjectionRef())

    ds_to_wgs84 = osr.CoordinateTransformation(
            ds_sr, wgs84_sr)

    gt = ds.GetGeoTransform()

    minx = gt[0]
    miny = gt[3] + ds.RasterXSize * gt[4] + ds.RasterYSize * gt[5]
    maxx = gt[0] + ds.RasterXSize * gt[1] + ds.RasterYSize * gt[2]
    maxy = gt[3]

    def t(x, y):
        lng, lat, _ = ds_to_wgs84.TransformPoint(x, y)
        return lat, lng

    return np.array([
        t(minx, miny),
        t(maxx, miny),
        t(maxx, maxy),
        t(minx, maxy)
    ]) 

def makeUtmCS(lon, lat):
    utmZone = utmGetZone(lon)
    isNorthern = utmIsNorthern(lat)
    # set utm coordinate system
    utmCs = osr.SpatialReference()
    utmCs.SetWellKnownGeogCS('WGS84')
    utmCs.SetUTM(utmZone,isNorthern)
    return utmCs

#########################################
# function to produce polygons from points
# inPoints is a list of lists e.g. [[[x1,y1], [x2,y2]], [[x3,y3], [x4,y4]]]
# each list of points is saved as a seperate feauture in the final file
######################################### 

def transPoint(point, inproj, outproj):
    sr_srs = osr.SpatialReference()
    sr_srs.ImportFromEPSG(inproj)
    dst_srs = osr.SpatialReference()
    dst_srs.ImportFromEPSG(outproj)

    coordTransform = osr.CoordinateTransformation(sr_srs, dst_srs)

    gps_point = ogr.Geometry(ogr.wkbPoint)
    gps_point.AddPoint(point[0], point[1])
    gps_point.Transform(coordTransform)
    x = float(gps_point.ExportToWkt().split()[1].split('(')[1])
    y = float(gps_point.ExportToWkt().split()[2])
    return [x,y] 

def localortho(lon, lat):
    """Create srs for local orthographic projection centered at lat, lon
    """
    local_srs = osr.SpatialReference()
    local_proj = '+proj=ortho +lat_0=%0.7f +lon_0=%0.7f +datum=WGS84 +units=m +no_defs ' % (lat, lon)
    local_srs.ImportFromProj4(local_proj)
    return local_srs

#Transform geometry to local orthographic projection, useful for width/height and area calc 

def get_ds_srs(ds):
    """Get srs object for GDAL Datset
    """
    ds_srs = osr.SpatialReference()
    ds_srs.ImportFromWkt(ds.GetProjectionRef())
    return ds_srs 

def __init__(self, bbox, epsg):
        self.bbox = bbox
        self.geom = bbox2geom(bbox)
        self.srs = osr.SpatialReference()
        self.srs.ImportFromEPSG(epsg)

#This provides a preference order for projections 

def convert_and_import(xml_file):
    QgsProject.instance().clear()
    with tempfile.NamedTemporaryFile(delete=True) as f:
        out_f = f.name
    config_file = os.path.join(os.path.dirname(__file__), "gmlasconf.xml")
    gdal.SetConfigOption("OGR_SQLITE_SYNCHRONOUS", "OFF")
    ds = gdal.OpenEx("GMLAS:{}".format(xml_file), open_options=['EXPOSE_METADATA_LAYERS=YES', 'CONFIG_FILE={}'.format(config_file)])
    srs = osr.SpatialReference()
    qgs_srs = QgsCoordinateReferenceSystem("EPSG:4326")
    srs.ImportFromWkt(qgs_srs.toWkt())
    params = {
        'destNameOrDestDS': out_f
        , 'srcDS': ds
        , 'format': "SQLite"
        , 'accessMode': "overwrite"
        , 'datasetCreationOptions': ['SPATIALITE=YES']
        , 'options' : ['-forceNullable', '-skipfailures']
        #, 'srcSRS': srs
        #, 'dstSRS': srs
        , 'geometryType': 'CONVERT_TO_LINEAR'
        , 'reproject': False
    }
    # call gdal to convert
    gdal.VectorTranslate(**params)
    # fix geometry types
    ds = None
    # populate the qgis project
    import_in_qgis(out_f, "SQLite")

    layers = []
    for lid in sorted(QgsProject.instance().mapLayers().keys()):
        vl = QgsProject.instance().mapLayer(lid)
        layers.append((vl.name(), vl.wkbType()))
    rels = []
    relations = QgsProject.instance().relationManager().relations()
    for relid in sorted(relations.keys()):
        rel = relations[relid]
        p = rel.fieldPairs()
        rels.append((rel.id()[0:3], rel.referencingLayer().name(), list(p.keys())[0], rel.referencedLayer().name(), list(p.values())[0]))

    return sorted(layers), sorted(rels) 

def dest_srs(self):
        srs = osr.SpatialReference()
        srs.ImportFromWkt(self.srsSelectionWidget.crs().toWkt())
        assert srs.Validate() == 0
        return srs 

def reproject_params(self, temp_datasource_path):
        params = {
            'destNameOrDestDS': temp_datasource_path,
            'srcDS': self.src_datasource(),
            'format': 'SQLite',
            'datasetCreationOptions': ['SPATIALITE=YES'],
            'dstSRS': self.dest_srs(),
            'reproject': True,
            'options': ['-skipfailures']
        }

        if self.bboxGroupBox.isChecked():
            if self.bboxWidget.value() == '':
                raise InputError("Extent is empty")
            if not self.bboxWidget.isValid():
                raise InputError("Extent is invalid")
            bbox = self.bboxWidget.rectangle()
            params['spatFilter'] = (bbox.xMinimum(),
                                    bbox.yMinimum(),
                                    bbox.xMaximum(),
                                    bbox.yMaximum())
            srs = osr.SpatialReference()
            srs.ImportFromWkt(self.bboxWidget.crs().toWkt())
            assert srs.Validate() == 0
            params['spatSRS'] = srs

        return params 

def src_srs(self):
        srs = osr.SpatialReference()
        srs.ImportFromWkt(self.sourceSrs.crs().toWkt())
        assert srs.Validate() == 0
        return srs 

def get_default_projection():
    """Create a default projection object (wgs84)"""
    proj = osr.SpatialReference()
    proj.ImportFromEPSG(4326)
    return proj 

def test_get_clip_mask(self):
        proj_gk = osr.SpatialReference()
        proj_gk.ImportFromEPSG(31466)
        coords = np.array([[2600020., 5630020.], [2600030., 5630030.],
                           [2600040., 5630040.], [2700100., 5630030.],
                           [2600040., 5640000.]])
        mask = zonalstats.get_clip_mask(coords, self.npobj, proj_gk)
        out = np.array([True, True, True, False, False])
        np.testing.assert_array_equal(mask, out) 

def test_proj4_to_osr(self):
        srs = georef.proj4_to_osr('+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 '
                                  '+k_0=0.99987742 +x_0=600000 +y_0=2200000 '
                                  '+a=6378249.2 +b=6356515 '
                                  '+towgs84=-168,-60,320,0,0,0,0 '
                                  '+pm=paris +units=m +no_defs')
        p4 = srs.ExportToProj4()
        srs2 = osr.SpatialReference()
        srs2.ImportFromProj4(p4)
        self.assertTrue(srs.IsSame(srs2)) 

def read_safnwc(filename):
    """Read MSG SAFNWC hdf5 file into a gdal georeferenced object

    Parameters
    ----------
    filename : string
        satellite file name

    Returns
    -------
    ds : gdal.DataSet
        with satellite data
    """

    root = gdal.Open(filename)
    ds1 = gdal.Open('HDF5:' + filename + '://CT')
    ds = gdal.GetDriverByName('MEM').CreateCopy('out', ds1, 0)

    # name = os.path.basename(filename)[7:11]
    try:
        proj = osr.SpatialReference()
        proj.ImportFromProj4(ds.GetMetadata()["PROJECTION"])
    except Exception:
        raise NameError("No metadata for satellite file %s" % filename)
    geotransform = root.GetMetadata()["GEOTRANSFORM_GDAL_TABLE"].split(",")
    geotransform[0] = root.GetMetadata()["XGEO_UP_LEFT"]
    geotransform[3] = root.GetMetadata()["YGEO_UP_LEFT"]
    ds.SetProjection(proj.ExportToWkt())
    ds.SetGeoTransform([float(x) for x in geotransform])
    return ds 

def convert_to_latlon(geoimg, lines):
    srs = osr.SpatialReference(geoimg.srs()).ExportToProj4()
    projin = Proj(srs)
    projout = Proj(init='epsg:4326')
    newlines = []
    for line in lines:
        l = []
        for point in line:
            pt = transform(projin, projout, point[0], point[1])
            l.append(pt)
        newlines.append(l)
    return antimeridian_linesplit(newlines) 

def __init__(self, filename):
        if os.path.isfile(filename):
            self.filename = filename if os.path.isabs(filename) else os.path.join(os.getcwd(), filename)
            self.raster = gdal.Open(filename, GA_ReadOnly)
        else:
            raise IOError('file does not exist')

        self.cols = self.raster.RasterXSize
        self.rows = self.raster.RasterYSize
        self.bands = self.raster.RasterCount
        self.dim = [self.rows, self.cols, self.bands]
        self.driver = self.raster.GetDriver()
        self.format = self.driver.ShortName
        self.dtype = gdal.GetDataTypeName(self.raster.GetRasterBand(1).DataType)
        self.projection = self.raster.GetProjection()
        self.srs = osr.SpatialReference(wkt=self.projection)
        self.proj4 = self.srs.ExportToProj4()
        try:
            self.epsg = crsConvert(self.srs, 'epsg')
        except RuntimeError:
            self.epsg = None
        self.geogcs = self.srs.GetAttrValue('geogcs')
        self.projcs = self.srs.GetAttrValue('projcs') if self.srs.IsProjected() else None
        self.geo = dict(
            zip(['xmin', 'xres', 'rotation_x', 'ymax', 'rotation_y', 'yres'], self.raster.GetGeoTransform()))

        # note: yres is negative!
        self.geo['xmax'] = self.geo['xmin'] + self.geo['xres'] * self.cols
        self.geo['ymin'] = self.geo['ymax'] + self.geo['yres'] * self.rows

        self.res = [abs(float(self.geo['xres'])), abs(float(self.geo['yres']))]
        self.nodata = self.raster.GetRasterBand(1).GetNoDataValue()

        self.__data = [] 

def crsConvert(crsIn, crsOut):
    """
    convert between epsg, wkt, proj4 and opengis spatial references
    crsText must be a osr.SpatialReference object, an opengis URL (e.g. 'http://www.opengis.net/def/crs/EPSG/0/4326') or a string of type WKT, PROJ4 or EPSG
    crsOut must be either 'wkt', 'proj4', 'epsg', 'osr', 'opengis' or 'prettyWkt' (a wkt string formatted for readability)
    if type 'osr' is selected the function will return a spatial reference object of type osr.SpatialReference()
    """
    if isinstance(crsIn, osr.SpatialReference):
        srs = crsIn.Clone()
    else:
        srs = osr.SpatialReference()
        try:
            if 'opengis.net/def/crs/EPSG/0/' in crsIn:
                crsIn = int(os.path.basename(crsIn.strip('/')))
            srs.ImportFromEPSG(crsIn)
        except (TypeError, RuntimeError):
            try:
                srs.ImportFromWkt(crsIn)
            except (TypeError, RuntimeError):
                try:
                    srs.ImportFromProj4(crsIn)
                except (TypeError, RuntimeError):
                    raise TypeError('crsText not recognized; must be of type WKT, PROJ4 or EPSG')
    if crsOut == 'wkt':
        return srs.ExportToWkt()
    elif crsOut == 'prettyWkt':
        return srs.ExportToPrettyWkt()
    elif crsOut == 'proj4':
        return srs.ExportToProj4()
    elif crsOut == 'epsg':
        srs.AutoIdentifyEPSG()
        return int(srs.GetAuthorityCode(None))
    elif crsOut == 'opengis':
        srs.AutoIdentifyEPSG()
        return 'http://www.opengis.net/def/crs/EPSG/0/{}'.format(srs.GetAuthorityCode(None))
    elif crsOut == 'osr':
        return srs
    else:
        raise ValueError('crsOut not recognized; must be either wkt, proj4, opengis or epsg') 

def convertNodesCoordinates(graph, sourceEPSG, targetEPSG):
	sourceProjection = osr.SpatialReference()
	sourceProjection.ImportFromEPSG(sourceEPSG)
	destinationProjection = osr.SpatialReference()
	destinationProjection.ImportFromEPSG(targetEPSG)
	coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	for node in graph.nodes():
		x, y = graph.node[node]['longitude'], graph.node[node]['latitude']
		point = ogr.Geometry(ogr.wkbPoint)
		point.AddPoint(x,y)
		point.Transform(coordTrans)
		graph.node[node]['longitude'], graph.node[node]['latitude'] = point.GetX(), point.GetY()
	return graph 

def convertCoordinates(coordinate, sourceEPSG, targetEPSG):
	sourceProjection = osr.SpatialReference()
	sourceProjection.ImportFromEPSG(sourceEPSG)
	destinationProjection = osr.SpatialReference()
	destinationProjection.ImportFromEPSG(targetEPSG)
	coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	point = ogr.Geometry(ogr.wkbPoint)
	lon, lat = coordinate
	point.AddPoint(lon,lat)
	point.Transform(coordTrans)
	return point.GetX(), point.GetY() 

def batchConvertCoordinates(sourceCoordinates, sourceEPSG, targetEPSG):
	targetCoordinates = []
	sourceProjection = osr.SpatialReference()
	sourceProjection.ImportFromEPSG(sourceEPSG)
	destinationProjection = osr.SpatialReference()
	destinationProjection.ImportFromEPSG(targetEPSG)
	coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	for lon, lat in sourceCoordinates:
		point = ogr.Geometry(ogr.wkbPoint)
		point.AddPoint(lon,lat)
		point.Transform(coordTrans)
		targetCoordinates.append((point.GetX(), point.GetY()))
	return targetCoordinates 

def toMgrs(latitude, longitude, precision=5):
    """ Converts geodetic (latitude and longitude) coordinates to an MGRS
    coordinate string, according to the current ellipsoid parameters.

    @param latitude - latitude value
    @param longitude - longitude value
    @param precision - precision level of MGRS string
    @returns - MGRS coordinate string
    """
    if math.fabs(latitude) > 90:
        raise MgrsException('Latitude outside of valid range (-90 to 90 degrees).')

    if (longitude < -180) or (longitude > 360):
        raise MgrsException('Longitude outside of valid range (-180 to 360 degrees).')

    if (precision < 0) or (precision > MAX_PRECISION):
        raise MgrsException('The precision must be between 0 and 5 inclusive.')

    hemisphere, zone, epsg = _epsgForWgs(latitude, longitude)
    src = osr.SpatialReference()
    src.ImportFromEPSG(4326)
    dst = osr.SpatialReference()
    dst.ImportFromEPSG(epsg)
    ct = osr.CoordinateTransformation(src, dst)
    x, y, z = ct.TransformPoint(longitude, latitude)

    if (latitude < -80) or (latitude > 84):
        # Convert to UPS
        mgrs = _upsToMgrs(hemisphere, x, y, precision)
    else:
        # Convert to UTM
        mgrs = _utmToMgrs(zone, hemisphere, latitude, longitude, x, y, precision)

    return mgrs 

def pixelToGeoCoord(xPix, yPix, inputRaster, sourceSR='', geomTransform='', targetSR=''):
    # If you want to garuntee lon lat output, specify TargetSR  otherwise, geocoords will be in image geo reference
    # targetSR = osr.SpatialReference()
    # targetSR.ImportFromEPSG(4326)
    # Transform can be performed at the polygon level instead of pixel level

    if targetSR =='':
        performReprojection=False
        targetSR = osr.SpatialReference()
        targetSR.ImportFromEPSG(4326)
    else:
        performReprojection=True

    if geomTransform=='':
        srcRaster = gdal.Open(inputRaster)
        geomTransform = srcRaster.GetGeoTransform()

        source_sr = osr.SpatialReference()
        source_sr.ImportFromWkt(srcRaster.GetProjectionRef())




    geom = ogr.Geometry(ogr.wkbPoint)
    xOrigin = geomTransform[0]
    yOrigin = geomTransform[3]
    pixelWidth = geomTransform[1]
    pixelHeight = geomTransform[5]

    xCoord = (xPix * pixelWidth) + xOrigin
    yCoord = (yPix * pixelHeight) + yOrigin
    geom.AddPoint(xCoord, yCoord)


    if performReprojection:
        if sourceSR=='':
            srcRaster = gdal.Open(inputRaster)
            sourceSR = osr.SpatialReference()
            sourceSR.ImportFromWkt(srcRaster.GetProjectionRef())
        coord_trans = osr.CoordinateTransformation(sourceSR, targetSR)
        geom.Transform(coord_trans)

    return (geom.GetX(), geom.GetY()) 

def __init__(self, fname):
        ds = ogr.Open(fname)
        assert ds is not None, "Couldn't open %s" % fname
        layer = ds.GetLayer()
        assert layer is not None

        WGS84_SR = osr.SpatialReference()
        WGS84_SR.ImportFromEPSG(4326)

        shape_sr = osr.SpatialReference()
        shape_sr.ImportFromWkt(layer.GetSpatialRef().ExportToWkt())

        transform = osr.CoordinateTransformation(shape_sr, WGS84_SR)

        # http://geoinformaticstutorial.blogspot.ch/2012/10/accessing-vertices-from-polygon-with.html
        polygons = {}
        for feature in layer:
            attr = feature.items()
            newattr = {}
            # some attributes contains unicode character. ASCIIfy everything
            # TODO: A bit brutal, but easy...
            for k, v in attr.items():
                newk = k.decode('ascii', errors='ignore')
                newattr[newk] = v
            attr = newattr

            geometry = feature.GetGeometryRef()
            assert geometry.GetGeometryName() == 'POLYGON'
            # A ring is a polygon in shapefiles
            ring = geometry.GetGeometryRef(0)
            assert ring.GetGeometryName() == 'LINEARRING'
            # The ring duplicates the last point, so for the polygon to be
            # closed, last point should equal first point
            npoints = ring.GetPointCount()
            points = [ring.GetPoint(i) for i in xrange(npoints)]
            points = [transform.TransformPoint(*p) for p in points]
            # third column is elevation - discard
            points = np.array(points)[:, :2]
            # swap (lng, lat) to (lat, lng)
            points = points[:, ::-1]
            assert np.allclose(points[-1], points[0])

            name = attr['name']
            polygons[name] = points

        self.polygons = polygons