Geographic Information Systems Asked on December 18, 2021
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()
I ended up using the python port of ogr2ogr that can be found here
Answered by notmenotme on December 18, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP