Python osgeo.osr.SpatialReference() 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 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 clip_raster_to_intersection(raster_to_clip_path, extent_raster_path, out_raster_path, is_landsat=False):
    """
    Clips one raster to the extent proivded by the other raster, and saves the result at out_raster_path.
    Assumes both raster_to_clip and extent_raster are in the same projection.
    Parameters
    ----------
    raster_to_clip_path
        The location of the raster to be clipped.
    extent_raster_path
        The location of the raster that will provide the extent to clip to
    out_raster_path
        A location for the finished raster
    """

    with TemporaryDirectory() as td:
        temp_aoi_path = os.path.join(td, "temp_clip.shp")
        get_extent_as_shp(extent_raster_path, temp_aoi_path)
        ext_ras = gdal.Open(extent_raster_path)
        proj = osr.SpatialReference(wkt=ext_ras.GetProjection())
        srs_id = int(proj.GetAttrValue('AUTHORITY', 1))
        clip_raster(raster_to_clip_path, temp_aoi_path, out_raster_path, srs_id, flip_x_y = is_landsat) 

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 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 calc_raster_terrain_fixed_elevation(crs, elevation, grid_size, raster_path, locator, x_max, x_min, y_max,
                                        y_min):
    # local variables:
    temp_shapefile = locator.get_temporary_file("terrain.shp")
    cols = int((x_max - x_min) / grid_size)
    rows = int((y_max - y_min) / grid_size)
    shapes = Polygon([[x_min, y_min], [x_max, y_min], [x_max, y_max], [x_min, y_max], [x_min, y_min]])
    geodataframe = Gdf(index=[0], crs=crs, geometry=[shapes])
    geodataframe.to_file(temp_shapefile)
    # 1) opening the shapefile
    source_ds = ogr.Open(temp_shapefile)
    source_layer = source_ds.GetLayer()
    target_ds = gdal.GetDriverByName('GTiff').Create(raster_path, cols, rows, 1, gdal.GDT_Float32)  ##COMMENT 2
    target_ds.SetGeoTransform((x_min, grid_size, 0, y_max, 0, -grid_size))  ##COMMENT 3
    # 5) Adding a spatial reference ##COMMENT 4
    target_dsSRS = osr.SpatialReference()
    target_dsSRS.ImportFromProj4(crs)
    target_ds.SetProjection(target_dsSRS.ExportToWkt())
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(-9999)  ##COMMENT 5
    gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[elevation])  ##COMMENT 6
    target_ds = None  # closing the file 

def setup_input_srs(input_dataset, options):
    """
    Determines and returns the Input Spatial Reference System (SRS) as an osr object and as a
    WKT representation

    Uses in priority the one passed in the command line arguments. If None, tries to extract them
    from the input dataset
    """

    input_srs = None
    input_srs_wkt = None

    if options.s_srs:
        input_srs = osr.SpatialReference()
        input_srs.SetFromUserInput(options.s_srs)
        input_srs_wkt = input_srs.ExportToWkt()
    else:
        input_srs_wkt = input_dataset.GetProjection()
        if not input_srs_wkt and input_dataset.GetGCPCount() != 0:
            input_srs_wkt = input_dataset.GetGCPProjection()
        if input_srs_wkt:
            input_srs = osr.SpatialReference()
            input_srs.ImportFromWkt(input_srs_wkt)

    return input_srs, input_srs_wkt 

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 test_dump_raster(self):
        proj = osr.SpatialReference()
        proj.ImportFromEPSG(31466)
        filename = util.get_wradlib_data_file("shapefiles/agger/" "agger_merge.shp")
        test = zonalstats.DataSource(filename, srs=proj)
        test.dump_raster(
            tempfile.NamedTemporaryFile(mode="w+b").name,
            driver="netCDF",
            pixel_size=100.0,
        )
        test.dump_raster(
            tempfile.NamedTemporaryFile(mode="w+b").name,
            driver="netCDF",
            pixel_size=100.0,
            attr="FID",
        ) 

def get_radar_projection(sitecoords):
    """Get the native radar projection which is an azimuthal equidistant projection
    centered at the site using WGS84.

    Parameters
    ----------
    sitecoords : a sequence of two floats
        the WGS84 lon / lat coordinates of the radar location

    Returns
    -------
    proj : osr.SpatialReference
        projection definition

    """
    proj = osr.SpatialReference()
    proj.SetProjCS("Unknown Azimuthal Equidistant")
    proj.SetAE(sitecoords[1], sitecoords[0], 0, 0)

    return proj 

def get_extent(coords):
    """Get the extent of 2d coordinates

    Parameters
    ----------
    coords : :class:`numpy:numpy.ndarray`
        coordinates array with shape (...,(x,y))

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

    xmin = coords[..., 0].min()
    xmax = coords[..., 0].max()
    ymin = coords[..., 1].min()
    ymax = coords[..., 1].max()

    return xmin, xmax, ymin, ymax 

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

    Parameters
    ----------
    dataset : gdal.Dataset
        raster image with georeferencing

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

    Examples
    --------

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

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

def write_geotiff(tiff_path, data, upper_left_lon, upper_left_lat, step, AREA_OR_POINT='Point'):
    """
    写入geotiff。默认pixel is point，默认WGS84
    """
    rows, columns = data.shape
    bands = 1
    pixel_width = step
    pixel_height = -step
    half_step = step / 2
    upper_left_x = upper_left_lon - half_step
    upper_left_y = upper_left_lat + half_step
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(tiff_path, columns, rows, bands, gdal.GDT_Float32)
    dataset.SetMetadataItem('AREA_OR_POINT', AREA_OR_POINT)
    dataset.SetGeoTransform([upper_left_x, pixel_width, 0,
                             upper_left_y, 0, pixel_height])
    # 获取地理坐标系统信息，用于选取需要的地理坐标系统
    sr = osr.SpatialReference()
    sr.SetWellKnownGeogCS('WGS84')
    dataset.SetProjection(sr.ExportToWkt())  # 给新建图层赋予投影信息
    dataset.GetRasterBand(1).WriteArray(data)
    dataset.FlushCache()  # 将数据写入硬盘
    dataset = None 

def setup_input_srs(input_dataset, options):
    """
    Determines and returns the Input Spatial Reference System (SRS) as an osr object and as a
    WKT representation

    Uses in priority the one passed in the command line arguments. If None, tries to extract them
    from the input dataset
    """

    input_srs = None
    input_srs_wkt = None

    if options.s_srs:
        input_srs = osr.SpatialReference()
        input_srs.SetFromUserInput(options.s_srs)
        input_srs_wkt = input_srs.ExportToWkt()
    else:
        input_srs_wkt = input_dataset.GetProjection()
        if not input_srs_wkt and input_dataset.GetGCPCount() != 0:
            input_srs_wkt = input_dataset.GetGCPProjection()
        if input_srs_wkt:
            input_srs = osr.SpatialReference()
            input_srs.ImportFromWkt(input_srs_wkt)

    return input_srs, input_srs_wkt 

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 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 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 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 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 reproject_image(in_raster, out_raster_path, new_projection,  driver = "GTiff",  memory = 2e3, do_post_resample=True):
    """
    Creates a new, reprojected image from in_raster using the gdal.ReprojectImage function.

    Parameters
    ----------
    in_raster
        Either a gdal.Raster object or a path to a raster
    out_raster_path
        The path to the new output raster.
    new_projection
        The new projection in .wkt
    driver
        The format of the output raster.
    memory
        The amount of memory to give to the reprojection.
    do_post_resample
        If set to false, do not resample the image back to the original projection.

    Notes
    -----
    The GDAL reprojection routine changes the size of the pixels by a very small amount; for example, a 10m pixel image
    can become a 10.002m pixel resolution image. To stop alignment issues,
    by deafult this function resamples the images back to their original resolution

    """
    if type(new_projection) is int:
        proj = osr.SpatialReference()
        proj.ImportFromEPSG(new_projection)
        new_projection = proj.ExportToWkt()
    log = logging.getLogger(__name__)
    log.info("Reprojecting {} to {}".format(in_raster, new_projection))
    if type(in_raster) is str:
        in_raster = gdal.Open(in_raster)
    res = in_raster.GetGeoTransform()[1]
    gdal.Warp(out_raster_path, in_raster, dstSRS=new_projection, warpMemoryLimit=memory, format=driver)
    # After warping, image has irregular gt; resample back to previous pixel size
    # TODO: Make this an option
    if do_post_resample:
        resample_image_in_place(out_raster_path, res)
    return out_raster_path 

def preprocess_sen2_images(l2_dir, out_dir, l1_dir, cloud_threshold=60, buffer_size=0, epsg=None,
                           bands=("B02", "B03", "B04", "B08"), 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)
                resample_image_in_place(out_mask_path, out_resolution)
            else:
                log.info("Moving images to {}".format(out_dir))
                shutil.move(temp_path, out_path)
                shutil.move(mask_path, out_mask_path)
                resample_image_in_place(out_mask_path, out_resolution) 

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)


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