Python osgeo.osr.SpatialReference() Examples

The following are code examples for showing how to use osgeo.osr.SpatialReference(). They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don't like. You can also save this page to your account.

Example 1
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 7 votes vote down vote up
def convert_pixgwktList_to_wgs84wktList(inputRaster, wktPolygonPixList):
    ## returns # [[GeoWKT, PixWKT], ...]
    wgs84WKTList=[]
    if os.path.isfile(inputRaster):
        srcRaster = gdal.Open(inputRaster)
        targetSR = osr.SpatialReference()
        targetSR.ImportFromWkt(srcRaster.GetProjectionRef())
        geomTransform = srcRaster.GetGeoTransform()

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

        wgs84WKTList.extend(geom_wkt_list)

    # [[GeoWKT, PixWKT], ...]
    return wgs84WKTList 
Example 2
Project: HistoricalMap   Author: lennepkade   File: function_historical_map.py    (license) View Source Project 6 votes vote down vote up
def polygonize(self,rasterTemp,outShp):
            
            sourceRaster = gdal.Open(rasterTemp)
            band = sourceRaster.GetRasterBand(1)
            driver = ogr.GetDriverByName("ESRI Shapefile")
            # If shapefile already exist, delete it
            if os.path.exists(outShp):
                driver.DeleteDataSource(outShp)
                
            outDatasource = driver.CreateDataSource(outShp)            
            # get proj from raster            
            srs = osr.SpatialReference()
            srs.ImportFromWkt( sourceRaster.GetProjectionRef() )
            # create layer with proj
            outLayer = outDatasource.CreateLayer(outShp,srs)
            # Add class column (1,2...) to shapefile
      
            newField = ogr.FieldDefn('Class', ogr.OFTInteger)
            outLayer.CreateField(newField)
            
            gdal.Polygonize(band, None,outLayer, 0,[],callback=None)  
            
            outDatasource.Destroy()
            sourceRaster=None
            band=None
                    
            
            ioShpFile = ogr.Open(outShp, update = 1)
            
            
            lyr = ioShpFile.GetLayerByIndex(0)
            
            lyr.ResetReading()    

            for i in lyr:
                lyr.SetFeature(i)
            # if area is less than inMinSize or if it isn't forest, remove polygon 
                if i.GetField('Class')!=1:
                    lyr.DeleteFeature(i.GetFID())        
            ioShpFile.Destroy()
                        
            return outShp 
Example 3
Project: qgis-latlontools-plugin   Author: NationalSecurityAgency   File: mgrs.py    (license) View Source Project 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
Project: chorospy   Author: spyrostheodoridis   File: rasterFunc.py    (license) View Source Project 6 votes vote down vote up
def raster2array(rasterfn):
    raster = gdal.Open(rasterfn)
    band = raster.GetRasterBand(1)
    nodata = band.GetNoDataValue()
    array = band.ReadAsArray()
    
    proj = raster.GetProjection()
    inproj = osr.SpatialReference()
    inproj.ImportFromWkt(proj)
    
    geoTransform = raster.GetGeoTransform()
    minx = geoTransform[0]
    maxy = geoTransform[3]
    maxx = minx + geoTransform[1]*raster.RasterXSize
    miny = maxy + geoTransform[5]*raster.RasterYSize
    extent =  [minx, maxx, miny, maxy]
    pixelSizeXY = [geoTransform[1], geoTransform[5]]
    del raster, band
    return [array, nodata, extent, inproj, pixelSizeXY]

#clip a raster by vector 
Example 5
Project: layman   Author: CCSS-CZ   File: __init__.py    (license) View Source Project 6 votes vote down vote up
def _get_raster_attributes(self,ds,attrs):
        """Collect raster attributes
        """

        attrs["type"] = "raster"

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

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

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

        return attrs 
Example 6
Project: wradlib   Author: wradlib   File: raster.py    (license) View Source Project 6 votes vote down vote up
def read_gdal_projection(dset):
    """Get a projection (OSR object) from a GDAL dataset.

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

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

    Examples
    --------

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

    """
    wkt = dset.GetProjection()
    srs = osr.SpatialReference()
    srs.ImportFromWkt(wkt)
    # src = None
    return srs 
Example 7
Project: wradlib   Author: wradlib   File: projection.py    (license) View Source Project 6 votes vote down vote up
def proj4_to_osr(proj4str):
    """Transform a proj4 string to an osr spatial reference object

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

    Examples
    --------

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

    """
    proj = None
    if proj4str:
        proj = osr.SpatialReference()
        proj.ImportFromProj4(proj4str)
    else:
        proj = get_default_projection()
    return proj 
Example 8
Project: wradlib   Author: wradlib   File: projection.py    (license) View Source Project 6 votes vote down vote up
def epsg_to_osr(epsg=None):
    """Create osr spatial reference object from EPSG number

    .. versionadded:: 0.6.0

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

    Returns
    -------
    proj : osr.SpatialReference
        GDAL/OSR object defining projection
    """
    proj = None
    if epsg:
        proj = osr.SpatialReference()
        proj.ImportFromEPSG(epsg)
    else:
        proj = get_default_projection()
    return proj 
Example 9
Project: wradlib   Author: wradlib   File: projection.py    (license) View Source Project 6 votes vote down vote up
def wkt_to_osr(wkt=None):
    """Create osr spatial reference object from WKT string

    .. versionadded:: 0.6.0

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

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

    """
    proj = None
    if wkt:
        proj = osr.SpatialReference()
        proj.ImportFromWkt(wkt)
    else:
        proj = get_default_projection()
    return proj 
Example 10
Project: wradlib   Author: wradlib   File: test_georef.py    (license) View Source Project 6 votes vote down vote up
def test_get_radolan_grid_equality(self):
        # create radolan projection osr object
        scale = (1. + np.sin(np.radians(60.))) / (1. + np.sin(np.radians(90.)))
        dwd_string = ('+proj=stere +lat_0=90 +lat_ts=90 +lon_0=10 '
                      '+k={0:10.8f} +x_0=0 +y_0=0 +a=6370040 +b=6370040 '
                      '+to_meter=1000 +no_defs'.format(scale))
        proj_stereo = georef.proj4_to_osr(dwd_string)

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

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

        # check source and target arrays for equality
        self.assertTrue(np.allclose(radolan_grid_ll, self.radolan_grid_ll))
        self.assertTrue(np.allclose(radolan_grid_xy, self.radolan_grid_xy)) 
Example 11
Project: policosm   Author: ComplexCity   File: addMetricDistanceToEdges.py    (license) View Source Project 6 votes vote down vote up
def addMetricDistanceToEdge(x1,y1,x2,y2, epsgCode):
	# we assume initial epsg is wsg84 (merctor projection)
	
	if epsgCode != 4326:
		sourceProjection = osr.SpatialReference()
		sourceProjection.ImportFromEPSG(4326)
		destinationProjection = osr.SpatialReference()
		destinationProjection.ImportFromEPSG(epsgCode) # https://epsg.io/2154
		coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	
	line = ogr.Geometry(ogr.wkbLineString)
	line.AddPoint(x1,y1)
	line.AddPoint(x2,y2)
	if epsgCode != 4326:
		line.Transform(coordTrans)
	length = line.Length()

	return length 
Example 12
Project: policosm   Author: ComplexCity   File: addMetricDistanceToEdges.py    (license) View Source Project 6 votes vote down vote up
def addMetricDistanceToEdges(graph, epsgCode):
	# we assume initial epsg is wsg84 (merctor projection)
	metricDistance = {}
	
	if epsgCode != 4326:
		sourceProjection = osr.SpatialReference()
		sourceProjection.ImportFromEPSG(4326)
		destinationProjection = osr.SpatialReference()
		destinationProjection.ImportFromEPSG(epsgCode) # https://epsg.io/2154
		coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	
	for edge in graph.edges():
		node1, node2 =  edge
		line = ogr.Geometry(ogr.wkbLineString)
		line.AddPoint(graph.node[node1]['longitude'], graph.node[node1]['latitude'])
		line.AddPoint(graph.node[node2]['longitude'], graph.node[node2]['latitude'])
		if epsgCode != 4326:
			line.Transform(coordTrans)
		length = line.Length()
		metricDistance[edge] = length
	nx.set_edge_attributes(graph, 'length', metricDistance)

	return graph 
Example 13
Project: policosm   Author: ComplexCity   File: testAddMetricDistanceToEdges.py    (license) View Source Project 6 votes vote down vote up
def addMetricDistanceToEdges(graph, epsgCode):
	# we assume initial epsg is wsg84 (merctor projection)
	metricDistance = {}
	
	if epsgCode != 4326:
		sourceProjection = osr.SpatialReference()
		sourceProjection.ImportFromEPSG(4326)
		destinationProjection = osr.SpatialReference()
		destinationProjection.ImportFromEPSG(epsgCode) # https://epsg.io/2154
		coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	
	for edge in graph.edges():
		node1, node2 =  edge
		line = ogr.Geometry(ogr.wkbLineString)
		line.AddPoint(graph.node[node1]['longitude'], graph.node[node1]['latitude'])
		line.AddPoint(graph.node[node2]['longitude'], graph.node[node2]['latitude'])
		if epsgCode != 4326:
			line.Transform(coordTrans)
		length = line.Length()
		metricDistance[edge] = length
	nx.set_edge_attributes(graph, 'length', metricDistance)

	return graph 
Example 14
Project: mgrspy   Author: boundlessgeo   File: mgrs.py    (license) View Source Project 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 15
Project: DRIP-SLIP   Author: NASA-DEVELOP   File: SLIP.py    (license) View Source Project 6 votes vote down vote up
def array2raster(newRasterFilename,rasterOrigin,pixelWidth,pixelHeight,array,dataType):
	cols=array.shape[1]
	rows=array.shape[0]
	originX=rasterOrigin[0]
	originY=rasterOrigin[1]
	driver=gdal.GetDriverByName('GTiff')
	outRaster = driver.Create(newRasterFilename, cols, rows, 1, dataType)
	outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
	outband = outRaster.GetRasterBand(1)
	outband.WriteArray(array)
	outRasterSRS = osr.SpatialReference()
	outRasterSRS.ImportFromEPSG(32645)# this is the EPSG code for Nepal, should be changed for other locations
	outRaster.SetProjection(outRasterSRS.ExportToWkt())
	outband.FlushCache()

#takes a series of rasters and clips them to minExtent, then returns them as numpy arrays 
Example 16
Project: unmixing   Author: arthur-e   File: lsma.py    (license) View Source Project 6 votes vote down vote up
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() 
Example 17
Project: lopocs   Author: Oslandia   File: database.py    (GNU Lesser General Public License v2.1) View Source Project 5 votes vote down vote up
def srs(self):
        sr = SpatialReference()
        sr.ImportFromEPSG(self.srsid)
        return sr.ExportToWkt() 
Example 18
Project: pfb-network-connectivity   Author: azavea   File: detect_utm_zone.py    (license) View Source Project 5 votes vote down vote up
def get_bbox(shapefile):
    """ Gets the bounding box of a shapefile in EPSG 4326.
        If shapefile is not in WGS84, bounds are reprojected.
    """
    driver = ogr.GetDriverByName('ESRI Shapefile')
    # 0 means read, 1 means write
    data_source = driver.Open(shapefile, 0)
    # Check to see if shapefile is found.
    if data_source is None:
        failure('Could not open %s' % (shapefile))
    else:
        layer = data_source.GetLayer()
        shape_bbox = layer.GetExtent()

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

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

        return shape_bbox 
Example 19
Project: DeepOSM   Author: trailbehind   File: geo_util.py    (license) View Source Project 5 votes vote down vote up
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)) 
Example 20
Project: DeepOSM   Author: trailbehind   File: geo_util.py    (license) View Source Project 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 21
Project: DeepOSM   Author: trailbehind   File: geo_util.py    (license) View Source Project 5 votes vote down vote up
def pixel_to_lat_lon_web_mercator(raster_dataset, col, row):
    """Convert a pixel on the raster_dataset to web mercator (epsg:3857)."""
    spatial_reference = osr.SpatialReference()
    spatial_reference.ImportFromWkt(raster_dataset.GetProjection())
    ds_spatial_reference_proj_string = spatial_reference.ExportToProj4()
    in_proj = Proj(ds_spatial_reference_proj_string)
    out_proj = Proj(init='epsg:3857')

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

    x2, y2 = transform(in_proj, out_proj, ulon, ulat)
    x2, y2 = out_proj(x2, y2, inverse=True)
    return x2, y2 
Example 22
Project: utilities   Author: SpaceNetChallenge   File: externalVectorProcessing.py    (license) View Source Project 5 votes vote down vote up
def buildTindex(rasterFolder, rasterExtention='.tif'):
    rasterList = glob.glob(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))
    print(rasterList)

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

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

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

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



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

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

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


    return gTindex, gTindexLayer 
Example 23
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 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 24
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 5 votes vote down vote up
def returnBoundBox(xOff, yOff, pixDim, inputRaster, targetSR='', pixelSpace=False):
    # Returns Polygon for a specific square defined by a center Pixel and
    # number of pixels in each dimension.
    # Leave targetSR as empty string '' or specify it as a osr.SpatialReference()
    # targetSR = osr.SpatialReference()
    # targetSR.ImportFromEPSG(4326)
    if targetSR == '':
        targetSR = osr.SpatialReference()
        targetSR.ImportFromEPSG(4326)
    xCord = [xOff - pixDim / 2, xOff - pixDim / 2, xOff + pixDim / 2,
             xOff + pixDim / 2, xOff - pixDim / 2]
    yCord = [yOff - pixDim / 2, yOff + pixDim / 2, yOff + pixDim / 2,
             yOff - pixDim / 2, yOff - pixDim / 2]

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

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

    poly.AddGeometry(ring)

    return poly 
Example 25
Project: PostgreSQL-GeoPackage   Author: EOX-A   File: gpkg-pg_dump.py    (license) View Source Project 5 votes vote down vote up
def create_gpkg(
    gpkg_name, proj_string, size=(1, 1), geotransform=[0, 1, 0, 0, 0, -1],
    creation_options=None
):
    if os.path.exists("%s.gpkg" % gpkg_name):
        sys.stderr.write(
            "ERROR: SQLite GeoPackage '%s.gpkg' already exists.\n" % gpkg_name
        )
        sys.exit(1)

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

        proj = osr.SpatialReference()
        res = proj.SetWellKnownGeogCS(proj_string)
        if res != 0:
            if proj_string[0:4] == 'EPSG':
                proj.ImportFromEPSG(int(proj_string[5:]))
        gpkg.SetProjection(proj.ExportToWkt())
        gpkg.SetGeoTransform(geotransform)
        gpkg = None
    except Exception as e:
        sys.stderr.write(
            "ERROR: Cannot create SQLite GeoPackage '%s.gpkg'. "
            "Error message was: '%s'.\n" % (gpkg_name, e.message)
        )
        sys.exit(1) 
Example 26
Project: qgis-latlontools-plugin   Author: NationalSecurityAgency   File: mgrs.py    (license) View Source Project 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
    """
    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 27
Project: rastercube   Author: terrai   File: gdal_utils.py    (license) View Source Project 5 votes vote down vote up
def latlng_bounding_box_from_ds(ds):
    """
    Given a GDAL dataset, computes its bounding box int lat/lng
    """
    wgs84_sr = osr.SpatialReference()
    wgs84_sr.ImportFromEPSG(4326)

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

    ds_to_wgs84 = osr.CoordinateTransformation(ds_sr, wgs84_sr)

    geot = ds.GetGeoTransform()

    def _xy2latlng(x, y):
        x_geo = geot[0] + x * geot[1] + y * geot[2]
        y_geo = geot[3] + x * geot[4] + y * geot[5]
        lng, lat, _ = ds_to_wgs84.TransformPoint(x_geo, y_geo)
        return lat, lng
    poly = [
        _xy2latlng(0, 0),
        _xy2latlng(ds.RasterXSize, 0),
        _xy2latlng(ds.RasterXSize, ds.RasterYSize),
        _xy2latlng(0, ds.RasterYSize)
    ]
    return poly 
Example 28
Project: rastercube   Author: terrai   File: gdal_utils.py    (license) View Source Project 5 votes vote down vote up
def latlng_bbox_from_ds(ds):
    wgs84_sr = osr.SpatialReference()
    wgs84_sr.ImportFromEPSG(4326)

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

    ds_to_wgs84 = osr.CoordinateTransformation(
            ds_sr, wgs84_sr)

    gt = ds.GetGeoTransform()

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

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

    return np.array([
        t(minx, miny),
        t(maxx, miny),
        t(maxx, maxy),
        t(minx, maxy)
    ]) 
Example 29
Project: chorospy   Author: spyrostheodoridis   File: vectorFunc.py    (license) View Source Project 5 votes vote down vote up
def makeUtmCS(lon, lat):
    utmZone = utmGetZone(lon)
    isNorthern = utmIsNorthern(lat)
    # set utm coordinate system
    utmCs = osr.SpatialReference()
    utmCs.SetWellKnownGeogCS('WGS84')
    utmCs.SetUTM(utmZone,isNorthern)
    return utmCs

#########################################
# function to produce polygons from points
# inPoints is a list of lists e.g. [[[x1,y1], [x2,y2]], [[x3,y3], [x4,y4]]]
# each list of points is saved as a seperate feauture in the final file
######################################### 
Example 30
Project: chorospy   Author: spyrostheodoridis   File: transFunc.py    (license) View Source Project 5 votes vote down vote up
def transPoint(point, inproj, outproj):
    sr_srs = osr.SpatialReference()
    sr_srs.ImportFromEPSG(inproj)
    dst_srs = osr.SpatialReference()
    dst_srs.ImportFromEPSG(outproj)

    coordTransform = osr.CoordinateTransformation(sr_srs, dst_srs)

    gps_point = ogr.Geometry(ogr.wkbPoint)
    gps_point.AddPoint(point[0], point[1])
    gps_point.Transform(coordTransform)
    x = float(gps_point.ExportToWkt().split()[1].split('(')[1])
    y = float(gps_point.ExportToWkt().split()[2])
    return [x,y] 
Example 31
Project: pygeotools   Author: dshean   File: geolib.py    (license) View Source Project 5 votes vote down vote up
def localortho(lon, lat):
    """Create srs for local orthographic projection centered at lat, lon
    """
    local_srs = osr.SpatialReference()
    local_proj = '+proj=ortho +lat_0=%0.7f +lon_0=%0.7f +datum=WGS84 +units=m +no_defs ' % (lat, lon)
    local_srs.ImportFromProj4(local_proj)
    return local_srs

#Transform geometry to local orthographic projection, useful for width/height and area calc 
Example 32
Project: pygeotools   Author: dshean   File: geolib.py    (license) View Source Project 5 votes vote down vote up
def get_ds_srs(ds):
    """Get srs object for GDAL Datset
    """
    ds_srs = osr.SpatialReference()
    ds_srs.ImportFromWkt(ds.GetProjectionRef())
    return ds_srs 
Example 33
Project: pygeotools   Author: dshean   File: geolib.py    (license) View Source Project 5 votes vote down vote up
def __init__(self, bbox, epsg):
        self.bbox = bbox
        self.geom = bbox2geom(bbox)
        self.srs = osr.SpatialReference()
        self.srs.ImportFromEPSG(epsg)

#This provides a preference order for projections 
Example 34
Project: gml_application_schema_toolbox   Author: BRGM   File: test_load_in_qgis.py    (license) View Source Project 5 votes vote down vote up
def convert_and_import(xml_file):
    QgsProject.instance().clear()
    with tempfile.NamedTemporaryFile(delete=True) as f:
        out_f = f.name
    config_file = os.path.join(os.path.dirname(__file__), "gmlasconf.xml")
    gdal.SetConfigOption("OGR_SQLITE_SYNCHRONOUS", "OFF")
    ds = gdal.OpenEx("GMLAS:{}".format(xml_file), open_options=['EXPOSE_METADATA_LAYERS=YES', 'CONFIG_FILE={}'.format(config_file)])
    srs = osr.SpatialReference()
    qgs_srs = QgsCoordinateReferenceSystem("EPSG:4326")
    srs.ImportFromWkt(qgs_srs.toWkt())
    params = {
        'destNameOrDestDS': out_f
        , 'srcDS': ds
        , 'format': "SQLite"
        , 'accessMode': "overwrite"
        , 'datasetCreationOptions': ['SPATIALITE=YES']
        , 'options' : ['-forceNullable', '-skipfailures']
        #, 'srcSRS': srs
        #, 'dstSRS': srs
        , 'geometryType': 'CONVERT_TO_LINEAR'
        , 'reproject': False
    }
    # call gdal to convert
    gdal.VectorTranslate(**params)
    # fix geometry types
    ds = None
    # populate the qgis project
    import_in_qgis(out_f, "SQLite")

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

    return sorted(layers), sorted(rels) 
Example 35
Project: gml_application_schema_toolbox   Author: BRGM   File: export_gmlas_panel.py    (license) View Source Project 5 votes vote down vote up
def dest_srs(self):
        srs = osr.SpatialReference()
        srs.ImportFromWkt(self.srsSelectionWidget.crs().toWkt())
        assert srs.Validate() == 0
        return srs 
Example 36
Project: gml_application_schema_toolbox   Author: BRGM   File: export_gmlas_panel.py    (license) View Source Project 5 votes vote down vote up
def reproject_params(self, temp_datasource_path):
        params = {
            'destNameOrDestDS': temp_datasource_path,
            'srcDS': self.src_datasource(),
            'format': 'SQLite',
            'datasetCreationOptions': ['SPATIALITE=YES'],
            'dstSRS': self.dest_srs(),
            'reproject': True,
            'options': ['-skipfailures']
        }

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

        return params 
Example 37
Project: gml_application_schema_toolbox   Author: BRGM   File: import_gmlas_panel.py    (license) View Source Project 5 votes vote down vote up
def src_srs(self):
        srs = osr.SpatialReference()
        srs.ImportFromWkt(self.sourceSrs.crs().toWkt())
        assert srs.Validate() == 0
        return srs 
Example 38
Project: wradlib   Author: wradlib   File: projection.py    (license) View Source Project 5 votes vote down vote up
def get_default_projection():
    """Create a default projection object (wgs84)"""
    proj = osr.SpatialReference()
    proj.ImportFromEPSG(4326)
    return proj 
Example 39
Project: wradlib   Author: wradlib   File: test_zonalstats.py    (license) View Source Project 5 votes vote down vote up
def test_get_clip_mask(self):
        proj_gk = osr.SpatialReference()
        proj_gk.ImportFromEPSG(31466)
        coords = np.array([[2600020., 5630020.], [2600030., 5630030.],
                           [2600040., 5630040.], [2700100., 5630030.],
                           [2600040., 5640000.]])
        mask = zonalstats.get_clip_mask(coords, self.npobj, proj_gk)
        out = np.array([True, True, True, False, False])
        np.testing.assert_array_equal(mask, out) 
Example 40
Project: wradlib   Author: wradlib   File: test_georef.py    (license) View Source Project 5 votes vote down vote up
def test_proj4_to_osr(self):
        srs = georef.proj4_to_osr('+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 '
                                  '+k_0=0.99987742 +x_0=600000 +y_0=2200000 '
                                  '+a=6378249.2 +b=6356515 '
                                  '+towgs84=-168,-60,320,0,0,0,0 '
                                  '+pm=paris +units=m +no_defs')
        p4 = srs.ExportToProj4()
        srs2 = osr.SpatialReference()
        srs2.ImportFromProj4(p4)
        self.assertTrue(srs.IsSame(srs2)) 
Example 41
Project: wradlib   Author: wradlib   File: gdal.py    (license) View Source Project 5 votes vote down vote up
def read_safnwc(filename):
    """Read MSG SAFNWC hdf5 file into a gdal georeferenced object

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

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

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

    # name = os.path.basename(filename)[7:11]
    try:
        proj = osr.SpatialReference()
        proj.ImportFromProj4(ds.GetMetadata()["PROJECTION"])
    except Exception:
        raise NameError("No metadata for satellite file %s" % filename)
    geotransform = root.GetMetadata()["GEOTRANSFORM_GDAL_TABLE"].split(",")
    geotransform[0] = root.GetMetadata()["XGEO_UP_LEFT"]
    geotransform[3] = root.GetMetadata()["YGEO_UP_LEFT"]
    ds.SetProjection(proj.ExportToWkt())
    ds.SetGeoTransform([float(x) for x in geotransform])
    return ds 
Example 42
Project: beachfront-py   Author: venicegeo   File: vectorize.py    (license) View Source Project 5 votes vote down vote up
def convert_to_latlon(geoimg, lines):
    srs = osr.SpatialReference(geoimg.srs()).ExportToProj4()
    projin = Proj(srs)
    projout = Proj(init='epsg:4326')
    newlines = []
    for line in lines:
        l = []
        for point in line:
            pt = transform(projin, projout, point[0], point[1])
            l.append(pt)
        newlines.append(l)
    return antimeridian_linesplit(newlines) 
Example 43
Project: pyroSAR   Author: johntruckenbrodt   File: raster.py    (license) View Source Project 5 votes vote down vote up
def __init__(self, filename):
        if os.path.isfile(filename):
            self.filename = filename if os.path.isabs(filename) else os.path.join(os.getcwd(), filename)
            self.raster = gdal.Open(filename, GA_ReadOnly)
        else:
            raise IOError('file does not exist')

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

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

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

        self.__data = [] 
Example 44
Project: pyroSAR   Author: johntruckenbrodt   File: auxil.py    (license) View Source Project 5 votes vote down vote up
def crsConvert(crsIn, crsOut):
    """
    convert between epsg, wkt, proj4 and opengis spatial references
    crsText must be a osr.SpatialReference object, an opengis URL (e.g. 'http://www.opengis.net/def/crs/EPSG/0/4326') or a string of type WKT, PROJ4 or EPSG
    crsOut must be either 'wkt', 'proj4', 'epsg', 'osr', 'opengis' or 'prettyWkt' (a wkt string formatted for readability)
    if type 'osr' is selected the function will return a spatial reference object of type osr.SpatialReference()
    """
    if isinstance(crsIn, osr.SpatialReference):
        srs = crsIn.Clone()
    else:
        srs = osr.SpatialReference()
        try:
            if 'opengis.net/def/crs/EPSG/0/' in crsIn:
                crsIn = int(os.path.basename(crsIn.strip('/')))
            srs.ImportFromEPSG(crsIn)
        except (TypeError, RuntimeError):
            try:
                srs.ImportFromWkt(crsIn)
            except (TypeError, RuntimeError):
                try:
                    srs.ImportFromProj4(crsIn)
                except (TypeError, RuntimeError):
                    raise TypeError('crsText not recognized; must be of type WKT, PROJ4 or EPSG')
    if crsOut == 'wkt':
        return srs.ExportToWkt()
    elif crsOut == 'prettyWkt':
        return srs.ExportToPrettyWkt()
    elif crsOut == 'proj4':
        return srs.ExportToProj4()
    elif crsOut == 'epsg':
        srs.AutoIdentifyEPSG()
        return int(srs.GetAuthorityCode(None))
    elif crsOut == 'opengis':
        srs.AutoIdentifyEPSG()
        return 'http://www.opengis.net/def/crs/EPSG/0/{}'.format(srs.GetAuthorityCode(None))
    elif crsOut == 'osr':
        return srs
    else:
        raise ValueError('crsOut not recognized; must be either wkt, proj4, opengis or epsg') 
Example 45
Project: policosm   Author: ComplexCity   File: convertNodesCoordinates.py    (license) View Source Project 5 votes vote down vote up
def convertNodesCoordinates(graph, sourceEPSG, targetEPSG):
	sourceProjection = osr.SpatialReference()
	sourceProjection.ImportFromEPSG(sourceEPSG)
	destinationProjection = osr.SpatialReference()
	destinationProjection.ImportFromEPSG(targetEPSG)
	coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	for node in graph.nodes():
		x, y = graph.node[node]['longitude'], graph.node[node]['latitude']
		point = ogr.Geometry(ogr.wkbPoint)
		point.AddPoint(x,y)
		point.Transform(coordTrans)
		graph.node[node]['longitude'], graph.node[node]['latitude'] = point.GetX(), point.GetY()
	return graph 
Example 46
Project: policosm   Author: ComplexCity   File: convertCoordinates.py    (license) View Source Project 5 votes vote down vote up
def convertCoordinates(coordinate, sourceEPSG, targetEPSG):
	sourceProjection = osr.SpatialReference()
	sourceProjection.ImportFromEPSG(sourceEPSG)
	destinationProjection = osr.SpatialReference()
	destinationProjection.ImportFromEPSG(targetEPSG)
	coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	point = ogr.Geometry(ogr.wkbPoint)
	lon, lat = coordinate
	point.AddPoint(lon,lat)
	point.Transform(coordTrans)
	return point.GetX(), point.GetY() 
Example 47
Project: policosm   Author: ComplexCity   File: convertCoordinates.py    (license) View Source Project 5 votes vote down vote up
def batchConvertCoordinates(sourceCoordinates, sourceEPSG, targetEPSG):
	targetCoordinates = []
	sourceProjection = osr.SpatialReference()
	sourceProjection.ImportFromEPSG(sourceEPSG)
	destinationProjection = osr.SpatialReference()
	destinationProjection.ImportFromEPSG(targetEPSG)
	coordTrans = osr.CoordinateTransformation(sourceProjection, destinationProjection)
	for lon, lat in sourceCoordinates:
		point = ogr.Geometry(ogr.wkbPoint)
		point.AddPoint(lon,lat)
		point.Transform(coordTrans)
		targetCoordinates.append((point.GetX(), point.GetY()))
	return targetCoordinates 
Example 48
Project: mgrspy   Author: boundlessgeo   File: mgrs.py    (license) View Source Project 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
    """
    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 49
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 4 votes vote down vote up
def pixelToGeoCoord(xPix, yPix, inputRaster, sourceSR='', geomTransform='', targetSR=''):
    # If you want to garuntee lon lat output, specify TargetSR  otherwise, geocoords will be in image geo reference
    # targetSR = osr.SpatialReference()
    # targetSR.ImportFromEPSG(4326)
    # Transform can be performed at the polygon level instead of pixel level

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

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

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




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

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


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

    return (geom.GetX(), geom.GetY()) 
Example 50
Project: rastercube   Author: terrai   File: regions.py    (license) View Source Project 4 votes vote down vote up
def __init__(self, fname):
        ds = ogr.Open(fname)
        assert ds is not None, "Couldn't open %s" % fname
        layer = ds.GetLayer()
        assert layer is not None

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

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

        transform = osr.CoordinateTransformation(shape_sr, WGS84_SR)

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

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

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

        self.polygons = polygons