2

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:

QGIS Project after setting CRS in Python

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:

QGIS Project after setting CRS manually

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?

6
  • @Kadir I have used the Project Properties Window and Layer Properties Window in the QGIS GUI to set the correct layer and project CRS. When I do this, the image in the Map View "deforms" to display the map correctly. I want to do the same thing using Python. I'm trying to use the setCrs() method to accomplish this, but they are not "deforming". Does that help clarify the issue? Would it help if I posted the entire script? Thanks for your help. – sajonainnominado Jan 21 at 21:48
  • When I use projectObj.setCrs(projectCrs), map in the canvas deforms. What is your QGIS version? – Kadir Şahbaz Jan 21 at 23:03
  • @Kadir How do I use PyQGIS to change the CRS of a project or a layer so that it shows on the Map View? – sajonainnominado Jan 21 at 23:23
  • projectObj.setCrs(projectCrs) changes the CRS of a project. – Kadir Şahbaz Jan 22 at 10:07
  • Do you refresh the map canvas (programattically or using zoom or pan tool) after setting CRS? – Kadir Şahbaz Jan 22 at 11:05
0

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.

5
  • So, the SHP file used to make the layers comes from the NYC Housing Development Website. The SHP file has no set CRS. It comes with a ESRI PRJ file that describes the CRS. However, everything I have read suggests that QGIS cannot read PRJ files. That's why I'm trying to set the CRS. Also, I just want to know how to do it with the QGIS Python API. – sajonainnominado Jan 21 at 22:30
  • 1) Where did you read that? QGIS does read PRJ files. It might fail on some but that would be a different issue. 2) Have you tried my suggestion? 3) If it does not work, can you link to one of the Shapefiles you are using? – bugmenot123 Jan 22 at 8:35
  • @ bugmenot123 Sorry, I was thinking of this post, which indeed says that QGIS can read PRJ files, but might not be able to convert the reference to the ArcGIS CRS database to the QGIS CRS database. The SHP file, with DBF, PRJ and SHX files, can be downloaded here. I meant to link to it in my last reply. – sajonainnominado Jan 26 at 15:11
  • Thanks! You didn't add any detail on that in the initial question but seeing the data makes it clear that QGIS does not properly process that PRJ file indeed. – bugmenot123 Jan 26 at 17:39
  • I updated the answer – bugmenot123 Jan 26 at 17:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.