Python osgeo.osr.CoordinateTransformation() Examples

The following are 30 code examples of osgeo.osr.CoordinateTransformation(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module osgeo.osr , or try the search function .
Example #1
Source Project: cloudless   Author: BradNeuberg   File: download_planetlabs.py    License: Apache License 2.0 7 votes vote down vote up
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 
Example #2
Source Project: unmixing   Author: arthur-e   File: utils.py    License: MIT License 6 votes vote down vote up
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) 
Example #3
Source Project: QGISFMV   Author: All4Gis   File: mgrs.py    License: GNU General Public License v3.0 6 votes vote down vote up
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 
Example #4
Source Project: DsgTools   Author: dsgoficial   File: spatialiteDb.py    License: GNU General Public License v2.0 6 votes vote down vote up
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() 
Example #5
Source Project: DsgTools   Author: dsgoficial   File: shapefileDb.py    License: GNU General Public License v2.0 6 votes vote down vote up
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() 
Example #6
Source Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 5 votes vote down vote up
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 
Example #7
Source Project: apls   Author: CosmiQ   File: apls_utils.py    License: Apache License 2.0 5 votes vote down vote up
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)


############################################################################### 
Example #8
Source Project: apls   Author: CosmiQ   File: apls_utils.py    License: Apache License 2.0 5 votes vote down vote up
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)


############################################################################### 
Example #9
Source Project: bitcoin-arbitrage-trading-bot   Author: mammuth   File: jqvmap.py    License: MIT License 5 votes vote down vote up
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) 
Example #10
Source Project: qgis-latlontools-plugin   Author: NationalSecurityAgency   File: mgrs.py    License: GNU General Public License v2.0 5 votes vote down vote up
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 
Example #11
Source Project: gdal2tiles   Author: Luqqk   File: gdal2tiles.py    License: MIT License 5 votes vote down vote up
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 
Example #12
Source Project: QGISFMV   Author: All4Gis   File: mgrs.py    License: GNU General Public License v3.0 5 votes vote down vote up
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 
Example #13
Source Project: cube-in-a-box   Author: opendatacube   File: ls_public_bucket.py    License: MIT License 5 votes vote down vote up
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()} 
Example #14
Source Project: gdal2tiles   Author: tehamalab   File: gdal2tiles.py    License: MIT License 5 votes vote down vote up
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 
Example #15
Source Project: FloodMapsWorkshop   Author: vightel   File: hand_floodfill.py    License: Apache License 2.0 5 votes vote down vote up
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)
	# 
Example #16
Source Project: FloodMapsWorkshop   Author: vightel   File: srtm2_dem.py    License: Apache License 2.0 5 votes vote down vote up
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 
Example #17
Source Project: FloodMapsWorkshop   Author: vightel   File: dem.py    License: Apache License 2.0 5 votes vote down vote up
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 
Example #18
Source Project: FloodMapsWorkshop   Author: vightel   File: hand.py    License: Apache License 2.0 5 votes vote down vote up
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)
	# 
Example #19
Source Project: DeepOSM   Author: trailbehind   File: geo_util.py    License: MIT License 5 votes vote down vote up
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) 
Example #20
Source Project: neon   Author: NervanaSystems   File: spacenet_utils.py    License: Apache License 2.0 5 votes vote down vote up
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) 
Example #21
Source Project: Start-MAJA   Author: CNES   File: lib_mnt.py    License: GNU General Public License v3.0 5 votes vote down vote up
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" 
Example #22
Source Project: SMAC-M   Author: LarsSchy   File: rasterutils.py    License: MIT License 5 votes vote down vote up
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 
Example #23
Source Project: SMAC-M   Author: LarsSchy   File: rasterutils.py    License: MIT License 5 votes vote down vote up
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 
Example #24
Source Project: unmixing   Author: arthur-e   File: utils.py    License: MIT License 4 votes vote down vote up
def pixel_to_xy(pixel_pairs, gt=None, wkt=None, path=None, dd=False):
    '''
    Modified from code by Zachary Bears (zacharybears.com/using-python-to-
    translate-latlon-locations-to-pixels-on-a-geotiff/).
    This method translates given pixel locations into longitude-latitude
    locations on a given GeoTIFF. Arguments:
        pixel_pairs The pixel pairings to be translated in the form
                    [[x1, y1],[x2, y2]]
        gt          [Optional] A GDAL GeoTransform tuple
        wkt         [Optional] Projection information as Well-Known Text
        path        The file location of the GeoTIFF
        dd          True to use decimal degrees for longitude-latitude (False
                    is the default)

    NOTE: This method does not take into account pixel size and assumes a
            high enough image resolution for pixel size to be insignificant.
    '''
    assert path is not None or (gt is not None and wkt is not None), 'Function requires either a reference dataset or a geotransform and projection'

    pixel_pairs = map(list, pixel_pairs)
    srs = osr.SpatialReference() # Create a spatial ref. for dataset

    if path is not None:
        ds = gdal.Open(path) # Load the image dataset
        gt = ds.GetGeoTransform() # Get geotransform of the dataset

        if dd:
            srs.ImportFromWkt(ds.GetProjection()) # Set up coord. transform.

    else:
        srs.ImportFromWkt(wkt)

    # Will use decimal-degrees if so specified
    if dd:
        ct = osr.CoordinateTransformation(srs, srs.CloneGeogCS())

    # Go through all the point pairs and translate them to pixel pairings
    xy_pairs = []
    for point in pixel_pairs:
        # Translate the pixel pairs into untranslated points
        lon = point[0] * gt[1] + gt[0]
        lat = point[1] * gt[5] + gt[3]
        if dd:
            (lon, lat, holder) = ct.TransformPoint(lon, lat) # Points to space

        xy_pairs.append((lon, lat)) # Add the point to our return array

    return xy_pairs 
Example #25
Source Project: unmixing   Author: arthur-e   File: utils.py    License: MIT License 4 votes vote down vote up
def xy_to_pixel(xy_pairs, gt=None, wkt=None, path=None, dd=False):
    '''
    Modified from code by Zachary Bears (zacharybears.com/using-python-to-
    translate-latlon-locations-to-pixels-on-a-geotiff/).
    This method translates given longitude-latitude pairs into pixel
    locations on a given GeoTIFF. Arguments:
        xy_pairs    The X-Y coordinate pairs (e.g., decimal lat/lon pairs) to
                    be translated in the form [[lon1, lat1], [lon2, lat2]]
        gt          [Optional] A GDAL GeoTransform tuple
        wkt         [Optional] Projection information as Well-Known Text
        path        The file location of the GeoTIFF
        dd          True to use decimal degrees for longitude-latitude (False
                    is the default)

    NOTE: This method does not take into account pixel size and assumes a
    high enough image resolution for pixel size to be insignificant.
    Essentially, xy_pairs assumed to be pixel centroids but may be assigned
    to the nearest neighbor of the intended pixel. For this reason, for
    relevant applications, target pixels should be embedded in a matrix
    of similar pixels, e.g., a 90-by-90 meter area, for Landsat pixels.
    '''
    assert path is not None or (gt is not None and wkt is not None), \
        'Function requires either a reference dataset or a geotransform and projection'

    xy_pairs = map(list, xy_pairs)
    srs = osr.SpatialReference() # Create a spatial ref. for dataset

    if path is not None:
        ds = gdal.Open(path) # Load the image dataset
        gt = ds.GetGeoTransform() # Get geotransform of the dataset

        if dd:
            srs.ImportFromWkt(ds.GetProjection()) # Set up coord. transform.

    else:
        srs.ImportFromWkt(wkt) # Set up coord. transform.

    # Will use decimal-degrees if so specified
    if dd:
        ct = osr.CoordinateTransformation(srs.CloneGeogCS(), srs)

    # Go through all the point pairs and translate them to lng-lat pairs
    pixel_pairs = []
    for point in xy_pairs:
        if dd:
            # Change the point locations into the GeoTransform space
            (point[0], point[1], holder) = ct.TransformPoint(point[0], point[1])

        # Translate the x and y coordinates into pixel values
        x = (point[0] - gt[0]) / gt[1]
        y = (point[1] - gt[3]) / gt[5]
        pixel_pairs.append((int(x), int(y))) # Add point to our return array

    return pixel_pairs 
Example #26
Source Project: unmixing   Author: arthur-e   File: lsma.py    License: MIT License 4 votes vote down vote up
def point_to_pixel_geometry(
        points, source_epsg=None, target_epsg=None, pixel_side_length=30):
    '''
    Where points is a list of X,Y tuples and X and Y are coordinates in
    meters, returns a series of OGR Polygons where each Polygon is the
    pixel extent with a given point at its center. Assumes square pixels.
    Arguments:
        points      Sequence of X,Y numeric pairs or OGR Point geometries
        source_epsg The EPSG code of the source projection (Optional)
        target_epsg The EPSG code of the target projection (Optional)
        pixel_side_length   The length of one side of the intended pixel
    '''
    polys = []
    # Convert points to atomic X,Y pairs if necessary
    if isinstance(points[0], ogr.Geometry):
        points = [(p.GetX(), p.GetY()) for p in points]

    source_ref = target_ref = None
    if all((source_epsg, target_epsg)):
        source_ref = osr.SpatialReference()
        target_ref = osr.SpatialReference()
        source_ref.ImportFromEPSG(source_epsg)
        target_ref.ImportFromEPSG(target_epsg)
        transform = osr.CoordinateTransformation(source_ref, target_ref)

    for p in points:
        r = pixel_side_length / 2 # Half the pixel width
        ring = ogr.Geometry(ogr.wkbLinearRing) # Create a ring
        vertices = [
            (p[0] - r, p[1] + r), # Top-left
            (p[0] + r, p[1] + r), # Top-right, clockwise from here...
            (p[0] + r, p[1] - r),
            (p[0] - r, p[1] - r),
            (p[0] - r, p[1] + r)  # Add top-left again to close ring
        ]

        # Coordinate transformation
        if all((source_ref, target_ref)):
            vertices = [transform.TransformPoint(*v)[0:2] for v in vertices]

        for vert in vertices:
            ring.AddPoint(*vert)
        poly = ogr.Geometry(ogr.wkbPolygon)
        poly.AddGeometry(ring)
        polys.append(poly)

    return polys 
Example #27
Source Project: apls   Author: CosmiQ   File: apls_utils.py    License: Apache License 2.0 4 votes vote down vote up
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())


############################################################################### 
Example #28
Source Project: pyTSEB   Author: hectornieto   File: PyTSEB.py    License: GNU General Public License v3.0 4 votes vote down vote up
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 
Example #29
Source Project: geoio   Author: DigitalGlobe   File: base.py    License: MIT License 4 votes vote down vote up
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) 
Example #30
Source Project: autoRIFT   Author: leiyangleon   File: GeogridOptical.py    License: Apache License 2.0 4 votes vote down vote up
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])]