Python osr.SpatialReference() Examples

def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array, nodata, EPSG):
    """This function take a regular array to create a raster with it"""
    print("I am dealing with nodata values")
    array[np.isnan(array)] = nodata # Dealing with Nodata values
    print("I am writing the raster")
    cols = array.shape[1]
    rows = array.shape[0]

    originX = rasterOrigin[0]
    originY = rasterOrigin[1]

    driver = gdal.GetDriverByName('ENVI')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float64)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    #outband.SetNoDataValue(nodata)
    outband.WriteArray(array)
    #outband.SetNoDataValue(nodata)

    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(EPSG)

    outRaster.SetProjection(outRasterSRS.ExportToWkt())

    outband.FlushCache() 

def get_epsg(prjfile):
    '''Get the epsg code from a projection file of a shapefile

    Args:
        prjfile: a .prj file of a shapefile

    Returns:
        str: EPSG code

    '''

    prj_file = open(prjfile, 'r')
    prj_txt = prj_file.read()
    srs = osr.SpatialReference()
    srs.ImportFromESRI([prj_txt])
    srs.AutoIdentifyEPSG()
    # return EPSG code
    return srs.GetAuthorityCode(None) 

def getproj(self, withgrid=False, projformat='pyproj'):
        from PseudoNetCDF.coordutil import getproj4_from_cf_var
        gridmapping = self.variables['crs']
        proj4str = getproj4_from_cf_var(gridmapping, withgrid=withgrid)
        preserve_units = withgrid
        if projformat == 'proj4':
            return proj4str
        elif projformat == 'pyproj':
            import pyproj
            # pyproj adds +units=m, which is not right for latlon/lonlat
            if '+proj=lonlat' in proj4str or '+proj=latlon' in proj4str:
                preserve_units = True
            return pyproj.Proj(proj4str, preserve_units=preserve_units)
        elif projformat == 'wkt':
            import osr
            srs = osr.SpatialReference()
            # Imports WKT to Spatial Reference Object
            srs.ImportFromProj4(proj4str)
            srs.ExportToWkt()  # converts the WKT to an ESRI-compatible format
            return srs.ExportToWkt()
        else:
            raise ValueError('projformat must be pyproj, proj4 or wkt') 

def to_wkt(self, target_wkt):
        # If we're going from WGS84 -> Spherical Mercator, use PyProj, because
        #  there seems to be a bug in OGR that gives us an offset. (GDAL
        #  does fine, though.
        if target_wkt == self.wkt:
            return self
        import osr
        dstSpatialRef = osr.SpatialReference()
        dstSpatialRef.ImportFromEPSG(d_name_to_epsg[target_wkt])
        # dstSpatialRef.ImportFromWkt(d_name_to_wkt[target_wkt])
        srcSpatialRef = osr.SpatialReference()
        srcSpatialRef.ImportFromEPSG(d_name_to_epsg[self.wkt])
        # srcSpatialRef.ImportFromWkt(self.wkt_)
        coordTransform = osr.CoordinateTransformation(srcSpatialRef, dstSpatialRef)
        a, b, c = coordTransform.TransformPoint(self.lon, self.lat)
        return Point(b, a, wkt=target_wkt) 

def Save_as_MEM(data='', geo='', projection=''):
    """
    This function save the array as a memory file

    Keyword arguments:
    data -- [array], dataset of the geotiff
    geo -- [minimum lon, pixelsize, rotation, maximum lat, rotation,
            pixelsize], (geospatial dataset)
    projection -- interger, the EPSG code
    """
    # save as a geotiff
    driver = gdal.GetDriverByName("MEM")
    dst_ds = driver.Create('', int(data.shape[1]), int(data.shape[0]), 1,
                           gdal.GDT_Float32, ['COMPRESS=LZW'])
    srse = osr.SpatialReference()
    if projection == '':
        srse.SetWellKnownGeogCS("WGS84")
    else:
        srse.SetWellKnownGeogCS(projection)
    dst_ds.SetProjection(srse.ExportToWkt())
    dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
    dst_ds.SetGeoTransform(geo)
    dst_ds.GetRasterBand(1).WriteArray(data)
    return(dst_ds) 

def loadMetadata(self):
        # open the raster and its spatial reference
        self.src = gdal.Open(self.tif_path)

        if self.src is None:
            raise Exception('Could not load GDAL file "%s"' % self.tif_path)
        spatial_reference_raster = osr.SpatialReference(self.src.GetProjection())

        # get the WGS84 spatial reference
        spatial_reference = osr.SpatialReference()
        spatial_reference.ImportFromEPSG(4326)  # WGS84

        # coordinate transformation
        self.coordinate_transform = osr.CoordinateTransformation(spatial_reference, spatial_reference_raster)
        gt = self.geo_transform = self.src.GetGeoTransform()
        dev = (gt[1] * gt[5] - gt[2] * gt[4])
        self.geo_transform_inv = (gt[0], gt[5] / dev, -gt[2] / dev,
                                  gt[3], -gt[4] / dev, gt[1] / dev) 

def make_geotiff(data, length, width, latn_p, lonw_p, dlat, dlon, outfile, compress_option):

    if data.dtype == np.float32:
        dtype = gdal.GDT_Float32
        nodata = np.nan  ## or 0?
    elif data.dtype == np.uint8:
        dtype = gdal.GDT_Byte
        nodata = None

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(outfile, width, length, 1, dtype, options=compress_option)
    outRaster.SetGeoTransform((lonw_p, dlon, 0, latn_p, 0, dlat))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(data)
    if nodata is not None: outband.SetNoDataValue(nodata)
    outRaster.SetMetadataItem('AREA_OR_POINT', 'Point')
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()

    return


#%% Main 

def Spatial_Reference(epsg, return_string=True):
    """
    Obtain a spatial reference from the EPSG parameter
    """
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg)
    if return_string:
        return srs.ExportToWkt()
    else:
        return srs 

def transform_shape(shape, ssrs, tsrs):
    """ Transform shape from ssrs to tsrs (all wkt) and return as wkt """
    ogrgeom = CreateGeometryFromWkt(shape)
    trans = CoordinateTransformation(SpatialReference(ssrs), SpatialReference(tsrs))
    ogrgeom.Transform(trans)
    wkt = ogrgeom.ExportToWkt()
    ogrgeom = None
    return wkt 

def Save_as_tiff(name='', data='', geo='', projection=''):
    """
    This function save the array as a geotiff

    Keyword arguments:
    name -- string, directory name
    data -- [array], dataset of the geotiff
    geo -- [minimum lon, pixelsize, rotation, maximum lat, rotation,
            pixelsize], (geospatial dataset)
    projection -- integer, the EPSG code
    """
    # save as a geotiff
    driver = gdal.GetDriverByName("GTiff")
    dst_ds = driver.Create(name, int(data.shape[1]), int(data.shape[0]), 1,
                           gdal.GDT_Float32, ['COMPRESS=LZW'])
    srse = osr.SpatialReference()
    if projection == '':
        srse.SetWellKnownGeogCS("WGS84")

    else:
        try:
            if not srse.SetWellKnownGeogCS(projection) == 6:
                srse.SetWellKnownGeogCS(projection)
            else:
                try:
                    srse.ImportFromEPSG(int(projection))
                except:
                    srse.ImportFromWkt(projection)
        except:
            try:
                srse.ImportFromEPSG(int(projection))
            except:
                srse.ImportFromWkt(projection)

    dst_ds.SetProjection(srse.ExportToWkt())
    dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
    dst_ds.SetGeoTransform(geo)
    dst_ds.GetRasterBand(1).WriteArray(data)
    dst_ds = None
    return() 

def layer_srs(self):
        if self._layer_srs is None:
            try:
                raster = gdal.Open(self.filepath)
                srs = osr.SpatialReference()
                srs.ImportFromWkt(raster.GetProjectionRef())
                self._layer_srs = srs.ExportToProj4()
            except RuntimeError:
                self._layer_srs = \
                    """+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs """

        return self._layer_srs 

def export_raster(self, raster, fname, title, res=None):
        """ Write array to raster file with gdal

        Parameters
        ----------
        raster :    ndarray
                    raster to be exported
        fname :     str
                    file name
        title :     str
                    gdal title of the file
        res :       int/float, optional
                    resolution of the raster in m, if not provided the same as
                    the input CHM
        """
        res = res if res else self.resolution

        fname.parent.mkdir(parents=True, exist_ok=True)

        driver = gdal.GetDriverByName('GTIFF')
        y_pixels, x_pixels = raster.shape
        gdal_file = driver.Create(
            f'{fname}', x_pixels, y_pixels, 1, gdal.GDT_Float32
        )
        gdal_file.SetGeoTransform(
            (self.ul_lon, res, 0., self.ul_lat, 0., -res)
        )
        dataset_srs = gdal.osr.SpatialReference()
        dataset_srs.ImportFromEPSG(self.epsg)
        gdal_file.SetProjection(dataset_srs.ExportToWkt())
        band = gdal_file.GetRasterBand(1)
        band.SetDescription(title)
        band.SetNoDataValue(0.)
        band.WriteArray(raster)
        gdal_file.FlushCache()
        gdal_file = None 

def convertGeometry(geom, srs):
    import ogr
    import osr
    in_ref = osr.SpatialReference()
    in_ref.SetFromUserInput(srs)
    out_ref = osr.SpatialReference()
    out_ref.SetFromUserInput('EPSG:4326')

    g = ogr.CreateGeometryFromJson(json.dumps(geom))
    g.AssignSpatialReference(in_ref)
    g.TransformTo(out_ref)
    return json.loads(g.ExportToJson()) 

def transform_to_wgs(getLong, getLat, EPSGa):
    source = osr.SpatialReference()
    source.ImportFromEPSG(EPSGa)

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

    transform = osr.CoordinateTransformation(source, target)

    point = ogr.CreateGeometryFromWkt("POINT (" + str(getLong[0]) + " " + str(getLat[0]) + ")")
    point.Transform(transform)
    return [point.GetX(), point.GetY()] 

def reproject_geometry(geom, inproj4, out_epsg):
    '''Reproject a wkt geometry based on EPSG code

    Args:
        geom (ogr-geom): an ogr geom objecct
        inproj4 (str): a proj4 string
        out_epsg (str): the EPSG code to which the geometry should transformed

    Returns
        geom (ogr-geometry object): the transformed geometry

    '''

    geom = ogr.CreateGeometryFromWkt(geom)
    # input SpatialReference
    spatial_ref_in = osr.SpatialReference()
    spatial_ref_in.ImportFromProj4(inproj4)

    # output SpatialReference
    spatial_ref_out = osr.SpatialReference()
    spatial_ref_out.ImportFromEPSG(int(out_epsg))

    # create the CoordinateTransformation
    coord_transform = osr.CoordinateTransformation(spatial_ref_in,
                                                   spatial_ref_out)
    try:
        geom.Transform(coord_transform)
    except:
        print(' ERROR: Not able to transform the geometry')
        sys.exit()

    return geom 

def epsg_to_wkt_projection(epsg_code):
    
    spatial_ref = osr.SpatialReference()
    spatial_ref.ImportFromEPSG(epsg_code)  
            
    return spatial_ref.ExpotToWkt() 

def array2raster(newRasterfn,rasterRef,rasterArray):  
    try:
        source_ds=gdal.Open(rasterRef)
    except RuntimeError:
        print("Unable to open %s" % rasterRef,)
    print("..............................")
    print(rasterArray.shape)
    gt=source_ds.GetGeoTransform() #获取栅格的地理空间变换数据
    print(gt)    
    cols=rasterArray.shape[1] #获取列数量
    rows=rasterArray.shape[0] #获取行数量
    originX=gt[0] #获取起始点X值
    originY=gt[3] #获取起始点Y值
    pixelWidth=gt[1] #单元(cell)栅格宽
    pixelHeight=gt[5] #单元(cell)栅格高
    '''A:建立栅格'''
    driver=gdal.GetDriverByName('GTiff') #用于栅格输出的驱动实例化
    outRaster=driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float64) #建立输出栅格.Create(Driver self, char const * utf8_path, int xsize, int ysize, int bands=1, GDALDataType eType, char ** options=None) -> Dataset
    '''B:设置空间变换'''
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight)) #设置输出栅格的空间变换参数，或者直接使用gt，保持与参考栅格相同设置
    '''C:给栅格波段赋值'''
    outband=outRaster.GetRasterBand(1) #获取输出栅格的一个输出波段
    outband.WriteArray(rasterArray) #将数组写入该波段
    '''D:设置栅格坐标系(投影)'''
    outRasterSRS=osr.SpatialReference() #空间参考实例化
    outRasterSRS.ImportFromWkt(source_ds.GetProjectionRef()) #设置空间参考为参考栅格的投影值
    outRaster.SetProjection(outRasterSRS.ExportToWkt()) #设置输出栅格的投影
    '''E:清空缓存与关闭栅格'''
    outband.FlushCache() #清空缓存
    source_ds=None #关闭栅格 

def get_crs(obj):
    """Get a coordinate reference system from any georegistered object."""
    if isinstance(obj, gpd.GeoDataFrame):
        return _check_crs(obj.crs)
    elif isinstance(obj, rasterio.DatasetReader):
        return _check_crs(obj.crs)
    elif isinstance(obj, gdal.Dataset):
        # rawr
        return _check_crs(int(osr.SpatialReference(wkt=obj.GetProjection()).GetAttrValue(
            'AUTHORITY', 1)))
    else:
        raise TypeError("solaris doesn't know how to extract a crs from an "
                        "object of type {}".format(type(obj))) 

def get_proj4(prjfile):
    '''Get the proj4 string from a projection file of a shapefile

    Args:
        prjfile: a .prj file of a shapefile

    Returns:
        str: PROJ4 code

    '''

    prj_file = open(prjfile, 'r')
    prj_string = prj_file.read()

    # Lambert error
    if '\"Lambert_Conformal_Conic\"' in prj_string:

        print(' ERROR: It seems you used an ESRI generated shapefile'
              ' with Lambert Conformal Conic projection. ')
        print(' This one is not compatible with Open Standard OGR/GDAL'
              ' tools used here. ')
        print(' Reproject your shapefile to a standard Lat/Long projection'
              ' and try again')
        exit(1)

    srs = osr.SpatialReference()
    srs.ImportFromESRI([prj_string])
    return srs.ExportToProj4() 

def getprojwkt(ifile, withgrid=False):
    import osr
    proj4str = getproj4(ifile, withgrid=withgrid)

    srs = osr.SpatialReference()
    # Imports WKT to Spatial Reference Object
    srs.ImportFromProj4(proj4str)
    srs.ExportToWkt()  # converts the WKT to an ESRI-compatible format
    return srs.ExportToWkt() 

def get_attrs(self):
        from osr import SpatialReference
        attrs = self.ds.GetMetadata()
        try:
            sp = SpatialReference(wkt=self.ds.GetProjection())
            proj4 = sp.ExportToProj4()
        except:
            warn('Could not identify projection')
        else:
            attrs['proj4'] = proj4
        return FrozenOrderedDict(attrs) 

def get_srid_via_osr(projection_data):
    """Attempt to retrive the SRID using osr's API."""

    srs = osr.SpatialReference()
    srs.ImportFromESRI([projection_data])
    srs.AutoIdentifyEPSG()

    code = srs.GetAuthorityCode(None)
    if code:
        return int(code) 

def Spatial_Reference(epsg, return_string=True):
    """
    Obtain a spatial reference from the EPSG parameter
    """
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg)
    if return_string:
        return srs.ExportToWkt()
    else:
        return srs 

def mk_geotiff_obj(raster, fn, bands=1, gdal_data_type=gdal.GDT_Float32,
                   lat=[46, 45], lon=[-73, -72]):
    """
    Creates a new geotiff file objects using the WGS84 coordinate system, saves
    it to disk, and returns a handle to the python file object and driver

    Parameters
    ------------
    raster : array
        Numpy array of the raster data to be added to the object
    fn : str
        Name of the geotiff file
    bands : int (optional)
        See :py:func:`gdal.GetDriverByName('Gtiff').Create
    gdal_data : gdal.GDT_<type>
        Gdal data type (see gdal.GDT_...)
    lat : list
        northern lat, southern lat
    lon : list
        [western lon, eastern lon]
    """
    NNi, NNj = raster.shape
    driver = gdal.GetDriverByName('GTiff')
    obj = driver.Create(fn, NNj, NNi, bands, gdal_data_type)
    pixel_height = -np.abs(lat[0] - lat[1]) / (NNi - 1.0)
    pixel_width = np.abs(lon[0] - lon[1]) / (NNj - 1.0)
    obj.SetGeoTransform([lon[0], pixel_width, 0, lat[0], 0, pixel_height])
    srs = osr.SpatialReference()
    srs.SetWellKnownGeogCS('WGS84')
    obj.SetProjection(srs.ExportToWkt())
    obj.GetRasterBand(1).WriteArray(raster)
    return obj, driver 

def raster_to_WSG_and_UTM(raster_path, lat, lon):

    raster = gdal.Open(raster_path)
    source_projection_wkt = raster.GetProjection()
    inSRS_converter = osr.SpatialReference()
    inSRS_converter.ImportFromProj4(get_projected_coordinate_system(lat, lon))
    target_projection_wkt = inSRS_converter.ExportToWkt()
    new_raster = gdal.AutoCreateWarpedVRT(raster, source_projection_wkt, target_projection_wkt,
                                          gdal.GRA_NearestNeighbour)
    return new_raster 

def reproject_dataset(dataset, output_file_name, wkt_from="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs", epsg_to=32643, pixel_spacing=1000, file_format='GTiff'):
    '''
    :param dataset: Modis subset name use gdal.GetSubdatasets()
    :param output_file_name: file location along with proper extension
    :param wkt_from: Modis wkt (default)
    :param epsg_to: integer default(32643)
    :param pixel_spacing: in metres
    :param file_format: default GTiff
    :return:
    '''
    wgs84 = osr.SpatialReference()
    wgs84.ImportFromEPSG(epsg_to)
    modis = osr.SpatialReference()
    modis.ImportFromProj4(wkt_from)
    tx = osr.CoordinateTransformation(modis, wgs84)
    g = gdal.Open(dataset)
    geo_t = g.GetGeoTransform()
    print geo_t
    x_size = g.RasterXSize
    y_size = g.RasterYSize
    (ulx, uly, ulz) = tx.TransformPoint(geo_t[0], geo_t[3])
    (lrx, lry, lrz) = tx.TransformPoint(geo_t[0] + (geo_t[1]*x_size), geo_t[3]+ (geo_t[5]*y_size))
    mem_drv = gdal.GetDriverByName(file_format)
    dest = mem_drv.Create(output_file_name, int((lrx-ulx)/pixel_spacing), int((uly - lry)/pixel_spacing), 1, gdal.GDT_Float32)
    new_geo = ([ulx, pixel_spacing, geo_t[2], uly, geo_t[4], -pixel_spacing])
    dest.SetGeoTransform(new_geo)
    dest.SetProjection(wgs84.ExportToWkt())
    gdal.ReprojectImage(g, dest, modis.ExportToWkt(), wgs84.ExportToWkt(), gdal.GRA_Bilinear)
    print "reprojected" 

def __init__(self, srs=4326):
        self.srs = osr.SpatialReference()
        self.srs.ImportFromEPSG(srs) 

def test_composite_across_projections_meters():
    os.chdir(os.path.dirname(os.path.abspath(__file__)))
    try:
        os.remove(r"test_outputs/composite_test.tif")
    except FileNotFoundError:
        pass
    try:
        shutil.rmtree(r"test_outputs/reprojected")
    except FileNotFoundError:
        pass
    os.mkdir(r"test_outputs/reprojected")
    epsg = 32736
    proj = osr.SpatialReference()
    proj.ImportFromEPSG(epsg)
    projection = proj.ExportToWkt() # Refactor this terrible nonsense later

    test_data = [r"test_data/S2A_MSIL2A_20180703T073611_N0206_R092_T36MZE_20180703T094637.tif",
                 r"test_data/S2B_MSIL2A_20180728T073609_N0206_R092_T37MBV_20180728T114325.tif"]
    pyeo.raster_manipulation.reproject_image(test_data[0], r"test_outputs/reprojected/0.tif", projection)
    pyeo.raster_manipulation.reproject_image(test_data[1], r"test_outputs/reprojected/1.tif", projection)
    pyeo.raster_manipulation.reproject_image(pyeo.filesystem_utilities.get_mask_path(test_data[0]), r"test_outputs/reprojected/0.msk", projection)
    pyeo.raster_manipulation.reproject_image(pyeo.filesystem_utilities.get_mask_path(test_data[1]), r"test_outputs/reprojected/1.msk", projection)

    out_file = r"test_outputs/composite_test.tif"
    pyeo.raster_manipulation.composite_images_with_mask([
       r"test_outputs/reprojected/0.tif",
        r"test_outputs/reprojected/1.tif"],
        out_file)
    image = gdal.Open("test_outputs/composite_test.tif")
    assert image
    image_array = image.GetVirtualMemArray()
    assert image_array.max() > 1 

def test_composite_across_projections():
    os.chdir(os.path.dirname(os.path.abspath(__file__)))
    try:
        os.remove(r"test_outputs/composite_test.tif")
    except FileNotFoundError:
        pass
    try:
        shutil.rmtree(r"test_outputs/reprojected")
    except FileNotFoundError:
        pass
    os.mkdir(r"test_outputs/reprojected")
    epsg = 4326
    proj = osr.SpatialReference()
    proj.ImportFromEPSG(epsg)
    projection = proj.ExportToWkt() # Refactor this terrible nonsense later
    test_data = [r"test_data/S2A_MSIL2A_20180703T073611_N0206_R092_T36MZE_20180703T094637.tif",
                 r"test_data/S2B_MSIL2A_20180728T073609_N0206_R092_T37MBV_20180728T114325.tif"]
    pyeo.raster_manipulation.reproject_image(test_data[0], r"test_outputs/reprojected/0.tif", projection)
    pyeo.raster_manipulation.reproject_image(test_data[1], r"test_outputs/reprojected/1.tif", projection)
    pyeo.raster_manipulation.reproject_image(pyeo.filesystem_utilities.get_mask_path(test_data[0]), r"test_outputs/reprojected/0.msk", projection)
    pyeo.raster_manipulation.reproject_image(pyeo.filesystem_utilities.get_mask_path(test_data[1]), r"test_outputs/reprojected/1.msk", projection)

    out_file = r"test_outputs/composite_test.tif"
    pyeo.raster_manipulation.composite_images_with_mask([
       r"test_outputs/reprojected/0.tif",
        r"test_outputs/reprojected/1.tif"],
        out_file)
    image = gdal.Open("test_outputs/composite_test.tif")
    assert image
    image_array = image.GetVirtualMemArray()
    assert image_array.max() > 10 

def _generate_latlon_transformer(raster):
    native_projection = osr.SpatialReference()
    native_projection.ImportFromWkt(raster.GetProjection())
    latlon_projection = osr.SpatialReference()
    latlon_projection.ImportFromEPSG(4326)
    geotransform = raster.GetGeoTransform()
    return osr.CoordinateTransformation(native_projection, latlon_projection), geotransform