Python osgeo.osr.CoordinateTransformation() Examples

def reproject(geom, from_epsg, to_epsg):
    """
    Reproject the given geometry from the given EPSG code to another
    """
    # Note: this is currently only accurate for the U.S.
    source = osr.SpatialReference()
    source.ImportFromEPSG(from_epsg)

    target = osr.SpatialReference()
    target.ImportFromEPSG(to_epsg)

    transform = osr.CoordinateTransformation(source, target)

    geom.Transform(transform)

    return geom 

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 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 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
    """

    # To avoid precision issues, which appear when using more than 6 decimal places
    latitude = round(latitude, 6)
    longitude = round(longitude, 6)

    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 insertFrame(self, scale, mi, inom, frame):
        self.checkAndOpenDb()
        srid = self.findEPSG()
        geoSrid = QgsCoordinateReferenceSystem(int(srid)).geographicCRSAuthId().split(':')[-1]
        ogr.UseExceptions()
        outputDS = self.buildOgrDatabase()
        outputLayer=outputDS.GetLayerByName('public_aux_moldura_a')
        newFeat=ogr.Feature(outputLayer.GetLayerDefn())
        auxGeom = ogr.CreateGeometryFromWkb(frame)
        #set geographic srid from frame
        geoSrs = ogr.osr.SpatialReference()
        geoSrs.ImportFromEPSG(int(geoSrid))
        auxGeom.AssignSpatialReference(geoSrs)
        #reproject geom
        outSpatialRef = outputLayer.GetSpatialRef()
        coordTrans = osr.CoordinateTransformation(geoSrs, outSpatialRef)
        auxGeom.Transform(coordTrans)
        newFeat.SetGeometry(auxGeom)
        newFeat.SetField('mi', mi)
        newFeat.SetField('inom', inom)
        newFeat.SetField('escala', str(scale))
        out=outputLayer.CreateFeature(newFeat)
        outputDS.Destroy() 

def insertFrame(self, scale, mi, inom, frame):
        self.checkAndOpenDb()
        srid = self.findEPSG()
        geoSrid = QgsCoordinateReferenceSystem(int(srid)).geographicCRSAuthId().split(':')[-1]
        ogr.UseExceptions()
        outputDS = self.buildOgrDatabase()
        outputLayer=outputDS.GetLayerByName(self.getFrameLayerName())
        newFeat=ogr.Feature(outputLayer.GetLayerDefn())
        auxGeom = ogr.CreateGeometryFromWkb(frame)
        #set geographic srid from frame
        geoSrs = ogr.osr.SpatialReference()
        geoSrs.ImportFromEPSG(int(geoSrid))
        auxGeom.AssignSpatialReference(geoSrs)
        #reproject geom
        outSpatialRef = outputLayer.GetSpatialRef()
        coordTrans = osr.CoordinateTransformation(geoSrs, outSpatialRef)
        auxGeom.Transform(coordTrans)
        newFeat.SetGeometry(auxGeom)
        newFeat.SetField('mi', mi)
        newFeat.SetField('inom', inom)
        newFeat.SetField('escala', str(scale))
        out=outputLayer.CreateFeature(newFeat)
        outputDS.Destroy() 

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 metersToLatLng(self,ds,X,Y):
		srs = osr.SpatialReference()
		srs.ImportFromWkt(ds.GetProjection())
		srsLatLong = srs.CloneGeogCS()
		ct = osr.CoordinateTransformation(srs,srsLatLong)
		return ct.TransformPoint(X,Y)

	#
	# Returns drainage direction data (HydroSHEDS)
	# 

def metersToLatLng(self,ds,X,Y):
		srs = osr.SpatialReference()
		srs.ImportFromWkt(ds.GetProjection())
		srsLatLong = srs.CloneGeogCS()
		ct = osr.CoordinateTransformation(srs,srsLatLong)
		return ct.TransformPoint(X,Y)
		
	# create matching osm water layer 

def metersToLatLng(self,ds,X,Y):
		srs = osr.SpatialReference()
		srs.ImportFromWkt(ds.GetProjection())
		srsLatLong = srs.CloneGeogCS()
		ct = osr.CoordinateTransformation(srs,srsLatLong)
		return ct.TransformPoint(X,Y)
		
	# create matching osm water layer 

def metersToLatLng(self,ds,X,Y):
		srs = osr.SpatialReference()
		srs.ImportFromWkt(ds.GetProjection())
		srsLatLong = srs.CloneGeogCS()
		ct = osr.CoordinateTransformation(srs,srsLatLong)
		return ct.TransformPoint(X,Y)

	#
	# Returns Height data (Conditioned DEM or VOID-filled from HydroSHEDS)
	# 

def get_tile_swne(tile_job_info, options):
    if options.profile == 'mercator':
        mercator = GlobalMercator(tileSize=options.tile_size)
        tile_swne = mercator.TileLatLonBounds
    elif options.profile == 'geodetic':
        geodetic = GlobalGeodetic(options.tmscompatible, tileSize=options.tile_size)
        tile_swne = geodetic.TileLatLonBounds
    elif options.profile == 'raster':
        srs4326 = osr.SpatialReference()
        srs4326.ImportFromEPSG(4326)
        if tile_job_info.kml and tile_job_info.in_srs_wkt:
            in_srs = osr.SpatialReference()
            in_srs.ImportFromWkt(tile_job_info.in_srs_wkt)
            ct = osr.CoordinateTransformation(in_srs, srs4326)

            def rastertileswne(x, y, z):
                pixelsizex = (2 ** (tile_job_info.tmaxz - z) * tile_job_info.out_geo_trans[1])
                west = tile_job_info.out_geo_trans[0] + x * tile_job_info.tilesize * pixelsizex
                east = west + tile_job_info.tilesize * pixelsizex
                south = tile_job_info.ominy + y * tile_job_info.tilesize * pixelsizex
                north = south + tile_job_info.tilesize * pixelsizex
                if not tile_job_info.is_epsg_4326:
                    # Transformation to EPSG:4326 (WGS84 datum)
                    west, south = ct.TransformPoint(west, south)[:2]
                    east, north = ct.TransformPoint(east, north)[:2]
                return south, west, north, east

            tile_swne = rastertileswne
        else:
            tile_swne = lambda x, y, z: (0, 0, 0, 0)   # noqa
    else:
        tile_swne = lambda x, y, z: (0, 0, 0, 0)   # noqa

    return tile_swne 

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 TestLand(lon, lat):
    latlon = osr.SpatialReference()
    latlon.ImportFromEPSG(4326)

    # create a point

    pt = ogr.Geometry(ogr.wkbPoint)
    pt.SetPoint_2D(0, lon, lat)

    # read shapefile
    shapefile = "land_polygons_osm/simplified_land_polygons.shp"
    driver = ogr.GetDriverByName("ESRI Shapefile")
    dataSource = driver.Open(shapefile, 0)
    layer = dataSource.GetLayer()
    targetProj = layer.GetSpatialRef()
    land = False

    # conversion to shapefile projection
    transform = osr.CoordinateTransformation(latlon, targetProj)
    pt.Transform(transform)

    # search point in layers
    for feature in layer:
        geom = feature.GetGeometryRef()
        if geom.Contains(pt):
            land = True
            break

    return land


##################################### Lecture de fichier de parametres "Mot_clé=Valeur" 

def get_coords(geo_ref_points, spatial_ref):
    t = osr.CoordinateTransformation(spatial_ref, spatial_ref.CloneGeogCS())

    def transform(p):
        # GDAL 3 reverses coordinate order, because... standards
        if LON_LAT_ORDER:
            # GDAL 2.0 order
            lon, lat, z = t.TransformPoint(p['x'], p['y'])
        else:
            # GDAL 3.0 order
            lat, lon, z = t.TransformPoint(p['x'], p['y'])
        return {'lon': lon, 'lat': lat}

    return {key: transform(p) for key, p in geo_ref_points.items()} 

def calculate_total_extent(data_path, list_of_files, to_srs):
    total_extent = [sys.maxsize, sys.maxsize, -
                    sys.maxsize - 1, -sys.maxsize - 1]
    for f in list_of_files:
        name = f[1]
        sub_dir = f[0] or ""
        source = gdal.Open(os.path.join(data_path, sub_dir, name))
        from_sys = osr.SpatialReference()
        from_sys.ImportFromWkt(source.GetProjection())
        gt = source.GetGeoTransform()
        width = source.RasterXSize
        height = source.RasterYSize

        thisMinX = gt[0]
        thisMinY = gt[3] + width*gt[4] + height*gt[5]
        thisMaxX = gt[0] + width*gt[1] + height*gt[2]
        thisMaxY = gt[3]

        to_sys = osr.SpatialReference()
        to_sys.ImportFromEPSG(to_srs)
        transformation = osr.CoordinateTransformation(from_sys, to_sys)
        minXY = transformation.TransformPoint(thisMinX, thisMinY)
        maxXY = transformation.TransformPoint(thisMaxX, thisMaxY)

        if minXY[0] < total_extent[0]:
            total_extent[0] = minXY[0]
        if minXY[1] < total_extent[1]:
            total_extent[1] = minXY[1]
        if maxXY[0] > total_extent[2]:
            total_extent[2] = maxXY[0]
        if maxXY[1] > total_extent[3]:
            total_extent[3] = maxXY[1]
    return total_extent

# This method is to be used by clients of this class 

def calculate_total_extent(data_path, list_of_files, to_srs):
    total_extent = [sys.maxsize, sys.maxsize, -
                    sys.maxsize - 1, -sys.maxsize - 1]
    for f in list_of_files:
        name = f[1]
        sub_dir = f[0] or ""
        source = gdal.Open(os.path.join(data_path, sub_dir, name))
        from_sys = osr.SpatialReference()
        from_sys.ImportFromWkt(source.GetProjection())
        gt = source.GetGeoTransform()
        width = source.RasterXSize
        height = source.RasterYSize

        thisMinX = gt[0]
        thisMinY = gt[3] + width*gt[4] + height*gt[5]
        thisMaxX = gt[0] + width*gt[1] + height*gt[2]
        thisMaxY = gt[3]

        to_sys = osr.SpatialReference()
        to_sys.ImportFromEPSG(to_srs)
        transformation = osr.CoordinateTransformation(from_sys, to_sys)
        minXY = transformation.TransformPoint(thisMinX, thisMinY)
        maxXY = transformation.TransformPoint(thisMaxX, thisMaxY)

        if minXY[0] < total_extent[0]:
            total_extent[0] = minXY[0]
        if minXY[1] < total_extent[1]:
            total_extent[1] = minXY[1]
        if maxXY[0] > total_extent[2]:
            total_extent[2] = maxXY[0]
        if maxXY[1] > total_extent[3]:
            total_extent[3] = maxXY[1]
    return total_extent

# This method is to be used by clients of this class 

def get_tile_swne(tile_job_info, options):
    if options.profile == 'mercator':
        mercator = GlobalMercator()
        tile_swne = mercator.TileLatLonBounds
    elif options.profile == 'geodetic':
        geodetic = GlobalGeodetic(options.tmscompatible)
        tile_swne = geodetic.TileLatLonBounds
    elif options.profile == 'raster':
        srs4326 = osr.SpatialReference()
        srs4326.ImportFromEPSG(4326)
        if tile_job_info.kml and tile_job_info.in_srs_wkt:
            in_srs = osr.SpatialReference()
            in_srs.ImportFromWkt(tile_job_info.in_srs_wkt)
            ct = osr.CoordinateTransformation(in_srs, srs4326)

            def rastertileswne(x, y, z):
                pixelsizex = (2 ** (tile_job_info.tmaxz - z) * tile_job_info.out_geo_trans[1])
                west = tile_job_info.out_geo_trans[0] + x * tile_job_info.tilesize * pixelsizex
                east = west + tile_job_info.tilesize * pixelsizex
                south = tile_job_info.ominy + y * tile_job_info.tilesize * pixelsizex
                north = south + tile_job_info.tilesize * pixelsizex
                if not tile_job_info.is_epsg_4326:
                    # Transformation to EPSG:4326 (WGS84 datum)
                    west, south = ct.TransformPoint(west, south)[:2]
                    east, north = ct.TransformPoint(east, north)[:2]
                return south, west, north, east

            tile_swne = rastertileswne
        else:
            tile_swne = lambda x, y, z: (0, 0, 0, 0)   # noqa
    else:
        tile_swne = lambda x, y, z: (0, 0, 0, 0)   # noqa

    return tile_swne 

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 _latlon2pixel(lat, lon, input_raster='', targetsr='', geom_transform=''):
    '''
    Convert latitude, longitude coords to pixexl coords.
    From spacenet geotools
    '''

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

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

    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 _wmp2pixel(x, y, input_raster='', targetsr='', geom_transform=''):
    '''
    Convert wmp coords to pixexl coords.
    '''

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

    geom = ogr.Geometry(ogr.wkbPoint)
    geom.AddPoint(x, y)

    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 intersect_rect(self, config, data_source):
    transform = osr.CoordinateTransformation( data_source.layer.GetSpatialRef(), data_source.spatialRef )
    point1 = transform.TransformPoint(config['rect'][0], config['rect'][1])
    point2 = transform.TransformPoint(config['rect'][2], config['rect'][3])
    rect = shapely.geometry.box(point1[0], point1[1], point2[0], point2[1])
    for geometry in data_source.geometries:
      geometry.geom = geometry.geom.intersection(rect) 

def _transform_osr(x1, y1, epsg_src, epsg_dst, polar=False):
    src = osr.SpatialReference()
    # Check if we are using osgeo.osr linked against PROJ 6+
    # If so, input axis ordering needs honored per projection, even though
    #   OAMS_TRADITIONAL_GIS_ORDER should fix it (doesn't seem to work for UPS)
    # See GDAL/OGR migration guide for 2.4 to 3.0
    # https://github.com/OSGeo/gdal/blob/master/gdal/MIGRATION_GUIDE.TXT and
    # https://trac.osgeo.org/gdal/wiki/rfc73_proj6_wkt2_srsbarn#Axisorderissues
    osr_proj6 = hasattr(src, 'SetAxisMappingStrategy')
    if not polar and osr_proj6:
        src.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    src.ImportFromEPSG(epsg_src)
    _log_proj_crs(src, proj_desc='src', espg=epsg_src)
    dst = osr.SpatialReference()
    if not polar and osr_proj6:
        dst.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    dst.ImportFromEPSG(epsg_dst)
    _log_proj_crs(dst, proj_desc='dst', espg=epsg_dst)
    ct = osr.CoordinateTransformation(src, dst)
    if polar and osr_proj6:
        # only supported with osgeo.osr v3.0.0+
        y2, x2, _ = ct.TransformPoint(y1, x1)
    else:
        x2, y2, _ = ct.TransformPoint(x1, y1)

    return x2, y2 

def get_data_from_vec_extent(self, vector=None, **kwargs):
        """This is a convenience method to find the extent of a vector and
        return the data from that extent.  kwargs can be anything accepted
        by get_data."""
        if vector is None:
            raise ValueError("Requires a vector to read.  The vector can be " \
                             "a string that describes a vector object or a " \
                             "path to a valid vector file.")

        if 'window' in kwargs.keys():
            raise ValueError("The window argument is not valid for this " \
                             "method. The vector file passed in defines " \
                             "the retrieval geometry.")

        if 'geom' in kwargs.keys():
            raise ValueError("The geom argument is not valid for this " \
                             "method. The vector file passed in defines " \
                             "the retrieval geometry.")

        if 'mask' in kwargs.keys():
            raise ValueError("A mask request is not valid for this method " \
                             "because it retrives data from the full extent " \
                             "of the vector.  You might want a rasterize " \
                             "method or iter_vector?")

        # ToDo Test for overlap of geom and image data?

        obj = ogr.Open(vector)
        lyr = obj.GetLayer(0)
        lyr_sr = lyr.GetSpatialRef()

        img_proj = self.meta.projection_string
        img_trans = self.meta.geo_transform
        img_sr = osr.SpatialReference()
        img_sr.ImportFromWkt(img_proj)

        coord_trans = osr.CoordinateTransformation(lyr_sr,img_sr)

        extent = self.ogr_extent_to_extent(lyr.GetExtent())

        window = self.extent_to_window(extent,coord_trans)
        [xoff, yoff, win_xsize, win_ysize] = window

        return self.get_data(window = window, **kwargs) 

def _get_subset(self, roi_shape, raster_proj_wkt, raster_geo_transform):

        # Find extent of ROI in roiShape projection
        roi = ogr.Open(roi_shape)
        roi_layer = roi.GetLayer()
        roi_extent = roi_layer.GetExtent()

        # Convert the extent to raster projection
        roi_proj = roi_layer.GetSpatialRef()
        raster_proj = osr.SpatialReference()
        raster_proj.ImportFromWkt(raster_proj_wkt)
        transform = osr.CoordinateTransformation(roi_proj, raster_proj)
        point_UL = ogr.CreateGeometryFromWkt("POINT ({} {})"
                                             .format(min(roi_extent[0], roi_extent[1]),
                                                     max(roi_extent[2], roi_extent[3])))
        point_UL.Transform(transform)
        point_UL = point_UL.GetPoint()
        point_LR = ogr.CreateGeometryFromWkt("POINT ({} {})"
                                             .format(max(roi_extent[0], roi_extent[1]),
                                                     min(roi_extent[2], roi_extent[3])))
        point_LR.Transform(transform)
        point_LR = point_LR.GetPoint()

        # Get pixel location of this extent
        ulX = raster_geo_transform[0]
        ulY = raster_geo_transform[3]
        pixel_size = raster_geo_transform[1]
        pixel_UL = [max(int(math.floor((ulY - point_UL[1]) / pixel_size)), 0),
                    max(int(math.floor((point_UL[0] - ulX) / pixel_size)), 0)]
        pixel_LR = [int(round((ulY - point_LR[1]) / pixel_size)),
                    int(round((point_LR[0] - ulX) / pixel_size))]

        # Get projected extent
        point_proj_UL = (ulX + pixel_UL[1] * pixel_size, ulY - pixel_UL[0] * pixel_size)

        # Convert to xoff, yoff, xcount, ycount as required by GDAL ReadAsArray()
        subset_pix = [pixel_UL[1], pixel_UL[0],
                      pixel_LR[1] - pixel_UL[1], pixel_LR[0] - pixel_UL[0]]

        # Get the geo transform of the subset
        subset_geo_transform = [point_proj_UL[0], pixel_size, raster_geo_transform[2],
                                point_proj_UL[1], raster_geo_transform[4], -pixel_size]

        return subset_pix, subset_geo_transform 

def determineBbox(self, zrange=[-200,4000]):
        '''
        Dummy.
        '''
        import numpy as np
        import datetime
        from osgeo import osr
        
#        import pdb
#        pdb.set_trace()


        samples = self.startingX + np.array([0, self.numberOfSamples]) * self.XSize
        lines = self.startingY + np.array([0, self.numberOfLines]) * self.YSize

        coordDat = osr.SpatialReference()
        if self.epsgDat:
            coordDat.ImportFromEPSG(self.epsgDat)
        else:
            raise Exception('EPSG code does not exist for image data')
    

        coordDem = osr.SpatialReference()
        if self.epsgDem:
            coordDem.ImportFromEPSG(self.epsgDem)
        else:
            raise Exception('EPSG code does not exist for DEM')
        
        
        trans = osr.CoordinateTransformation(coordDat, coordDem)
        
        

        utms = []
        xyzs = []


        ### Four corner coordinates
        for ss in samples:
            for ll in lines:
                for zz in zrange:
                    utms.append([ss,ll,zz])
                    x,y,z = trans.TransformPoint(ss, ll, zz)
                    xyzs.append([x,y,z])

        utms = np.array(utms)
        xyzs = np.array(xyzs)

        self._xlim = [np.min(xyzs[:,0]), np.max(xyzs[:,0])]
        self._ylim = [np.min(xyzs[:,1]), np.max(xyzs[:,1])] 

def translateLayer(self, inputLayer, inputLayerName, outputLayer, outputFileName, layerPanMap, errorDict, defaults={}, translateValues={}):
        """
        Makes the layer conversion
        """
        inputLayer.ResetReading()
        inSpatialRef = inputLayer.GetSpatialRef()
        outSpatialRef = outputLayer.GetSpatialRef()
        coordTrans = None
        if not inSpatialRef.IsSame(outSpatialRef):
            coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
        initialCount = outputLayer.GetFeatureCount()
        count = 0
        feat=inputLayer.GetNextFeature()
        #for feat in inputLayer:
        while feat:
            if not feat.geometry():
                continue
            inputId = feat.GetFID()
            if feat.geometry().GetGeometryCount() > 1:
                #Deaggregator
                for geom in feat.geometry():
                    newFeat=ogr.Feature(outputLayer.GetLayerDefn())
                    newFeat.SetFromWithMap(feat,True,layerPanMap)
                    auxGeom = ogr.Geometry(newFeat.geometry().GetGeometryType())
                    auxGeom.AssignSpatialReference(newFeat.geometry().GetSpatialReference())
                    auxGeom.AddGeometry(geom)
                    if coordTrans != None:
                        auxGeom.Transform(coordTrans)
                    newFeat.SetGeometry(auxGeom)
                    out=outputLayer.CreateFeature(newFeat)
                    if out != 0:
                        self.utils.buildNestedDict(errorDict, [inputLayerName], [inputId])
                    else:
                        count += 1
            else:
                newFeat=ogr.Feature(outputLayer.GetLayerDefn())
                newFeat.SetFromWithMap(feat,True,layerPanMap)
                if coordTrans != None:
                    geom = feat.GetGeometryRef()
                    geom.Transform(coordTrans)
                    newFeat.SetGeometry(geom)
                out=outputLayer.CreateFeature(newFeat)
                if out != 0:
                    self.utils.buildNestedDict(errorDict, [inputLayerName], [inputId])
                else:
                    count += 1
            feat=inputLayer.GetNextFeature()
        return count 

def vector2tiles(cls, vector, pcov=0.0, ptile=0.0, tilelist=None):
        """ Return matching tiles and coverage % for provided vector """
        from osgeo import ogr, osr

        # open tiles vector
        v = open_vector(cls.get_setting('tiles'))
        shp = ogr.Open(v.Filename())
        if v.LayerName() == '':
            layer = shp.GetLayer(0)
        else:
            layer = shp.GetLayer(v.LayerName())

        # create and warp site geometry
        ogrgeom = ogr.CreateGeometryFromWkt(vector.WKT())
        srs = osr.SpatialReference(vector.Projection())
        trans = osr.CoordinateTransformation(srs, layer.GetSpatialRef())
        ogrgeom.Transform(trans)
        # convert to shapely
        geom = loads(ogrgeom.ExportToWkt())

        # find overlapping tiles
        tiles = {}
        layer.SetSpatialFilter(ogrgeom)
        layer.ResetReading()
        feat = layer.GetNextFeature()
        while feat is not None:
            tgeom = loads(feat.GetGeometryRef().ExportToWkt())
            if tgeom.intersects(geom):
                area = geom.intersection(tgeom).area
                if area != 0:
                    tile = cls.feature2tile(feat)
                    tiles[tile] = (area / geom.area, area / tgeom.area)
            feat = layer.GetNextFeature()

        # remove any tiles not in tilelist or that do not meet thresholds for % cover
        remove_tiles = []
        if tilelist is None:
            tilelist = tiles.keys()
        for t in tiles:
            if (tiles[t][0] < (pcov / 100.0)) or (tiles[t][1] < (ptile / 100.0)) or t not in tilelist:
                remove_tiles.append(t)
        for t in remove_tiles:
            tiles.pop(t, None)
        return tiles 

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 pixelToGeoCoord(xPix, yPix, inputRaster, sourceSR='', geomTransform='',
                    targetSR=''):
    '''From spacenet geotools'''
    # If you want to gauruntee 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())


###############################################################################