Project vector in epsg 3347 to wkt in Python using GDAL

Question

I am new to working with geospatial data. I am currently using gdal and rasterio to work in Python. I have a raster with the following crs definition:

GEOGCS["Coordinate System imported from GRIB file",DATUM["unnamed",SPHEROID["Sphere",6371229,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]

I believe this is a WKT projection definition.

I also have a vector (geojson) in EPSG:3347 that I would like to project into the same crs as my raster.

I modified the below code from the GDAL cookbook, however it is not correctly projecting my vector. I would like to reproject the vector to the same crs as my raster and save the projected vector to my disk.

from osgeo import ogr, osr
import os
import gdal

driver = ogr.GetDriverByName('GeoJSON')

# get the input layer
inDataSet = driver.Open("B:/Canada-d8-100m/Canada-d8-100m.geojson")
inLayer = inDataSet.GetLayer()

inSpatialRef = inLayer.GetSpatialRef()

# loading projection
sr = osr.SpatialReference(str(inSpatialRef))

# detecting EPSG/SRID
res = sr.AutoIdentifyEPSG()

srid = sr.GetAuthorityCode(None)

print(srid)

# input SpatialReference
inSpatialRef.ImportFromEPSG(int(srid))

#read in wkt projection from raster
ds = gdal.Open("B:/dd.weather.gc.ca/model_gem_global/15km/grib2/lat_lon/00/000/CMC_glb_TMP_TGL_2_latlon.15x.15_2020052500_P000_Copy.grib2")
inSRS_wkt = ds.GetProjection()  
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromWkt(inSRS_wkt)

# create the CoordinateTransformation
coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

# create the output layer
outputShapefile = "B:/projected2.geojson"

if os.path.exists(outputShapefile):
    driver.DeleteDataSource(outputShapefile)

outDataSet = driver.CreateDataSource(outputShapefile)
outLayer = outDataSet.CreateLayer("ProjectedPoints2", geom_type=ogr.wkbPoint)

# add fields
inLayerDefn = inLayer.GetLayerDefn()

for i in range(0, inLayerDefn.GetFieldCount()):
    fieldDefn = inLayerDefn.GetFieldDefn(i)
    outLayer.CreateField(fieldDefn)

# get the output layer's feature definition
outLayerDefn = outLayer.GetLayerDefn()

# loop through the input features
inFeature = inLayer.GetNextFeature()

while inFeature:
    # get the input geometry
    geom = inFeature.GetGeometryRef()
    #print(geom)
    # reproject the geometry
    geom.Transform(coordTrans)
    #print(geom)
    # create a new feature
    outFeature = ogr.Feature(outLayerDefn)
    # set the geometry and attribute
    outFeature.SetGeometry(geom)
    for i in range(0, outLayerDefn.GetFieldCount()):
        outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
    # add the feature to the shapefile
    outLayer.CreateFeature(outFeature)
    # destroy the features and get the next input feature
    outFeature.Destroy()
    inFeature.Destroy()
    inFeature = inLayer.GetNextFeature()

# close the shapefiles
inDataSet.Destroy()
outDataSet.Destroy()

Answer 1

I ended up using the python port of ogr2ogr that can be found here