Setting project and layer CRS in PyQGIS with setCrs() has no effect on the map canvas

Question

I'm using the following code to set the project CRS in a Python script:
#creates the project and projectRoot objects
projectObj = QgsProject.instance()
projectRoot = projectObj.layerTreeRoot()

#sets the project CRS
projectCrs = QgsCoordinateReferenceSystem.fromEpsgId(3395)
projectObj.setCrs(projectCrs)
assert projectCrs.isValid()
print(projectCrs.isValid())

I'm using the following code in a for-loop to set the CRS for all layers:
shpSavedLayer.setCrs(projectCrs)

What I get is this:

Although, the QGIS Status Bar, Project Properties Window, and Layer Properties Window for every layer state that the CRS is set to EPSG:3395, there is clearly something wrong with the projection.
If I run the same code with the above lines commented out, then set the project CRS to EPSG:3395 via the Project Properties Window, I get this:

which looks a lot more accurate.
What's happening here? From the discussion on GIS SE, I'm guessing this has to do with the difference between setting the CRS and reprojecting the layer?

bugmenot123 · Answer

The line
shpSavedLayer.setCrs(projectCrs)

overrides your layers' CRSs regardless of their coordinates' actual CRS. For that to work properly you must know the actual CRS of your data.
The data you are using https://data.cityofnewyork.us/Housing-Development/Public-Use-Microdata-Areas-PUMA-/cwiz-gcty comes with a PRJ file of
GEOGCS[
  "WGS84(DD)", 
  DATUM["WGS84", SPHEROID["WGS84", 6378137.0, 298.257223563]],
  PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295], 
  AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH]
]

While this is actually a definition for WGS84 (EPSG:4326, at least the parameter values are the same), QGIS fails to properly process this as you must have noticed. But most importantly: That data is not in EPSG:3395.
To fix this, you need to set the layer's CRS separately:
layerCrs = QgsCoordinateReferenceSystem.fromEpsgId(4326)
shpSavedLayer.setCrs(layerCrs)

For the project/canvas you can use whatever CRS you like, EPSG:3395 8WGS 84 / World Mercator) would be a distorted view as well, I would recommend EPSG:26918 (NAD83 / UTM zone 18N) for optimal results.

Setting project and layer CRS in PyQGIS with setCrs() has no effect on the map canvas

One Answer

Add your own answers!

Ask a Question