Don't re-project projected polygon gdfs in metre-disaggregation by Evelyn-M · Pull Request #666

Don't re-project projected polygon gdfs in metre-disaggregation by Evelyn-M · Pull Request #666 · CLIMADA-project/climada_python

Changes proposed in this PR:

When disaggregating a polygon-containing geodataframe, one needs to specify whether resolution is specified in metres. For gdfs which are already in a projected crs (aka a metre-based one), it is currently strictly necessary to give the resolution (in metres), yet indicate to_metres=False (since else the code tries to re-project the already projected gdf, and fails). This is however a bit confusing to some end-users.
This hotfix checks in the mentioned case whether the gdf is already projected, and avoids re-projection.

This PR fixes #648

PR Author Checklist

Read the Contribution Guide
Correct target branch selected (if unsure, select develop)
Descriptive pull request title added
Source branch up-to-date with target branch
Documentation updated
Tests updated
Tests passing
No new linter issues
Changelog updated

PR Reviewer Checklist

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks a lot for a concise pull request and for working through the checklist! 😍 The change appears very sensible to me but I have trouble understanding the test. Maybe you can walk me through and we might find a way to simplify it a bit

Comment on lines +147 to +163

		#projected crs, to_meters=TRUE, FIX, dissag_val
		res = 20000
		EXP_POLY_PROJ = Exposures(GDF_POLY.to_crs(epsg=28992))
		exp_pnt = u_lp.exp_geom_to_pnt(
		EXP_POLY_PROJ, res=res, to_meters=True,
		disagg_met=u_lp.DisaggMethod.FIX, disagg_val=res**2
		)
		self.check_unchanged_exp(EXP_POLY_PROJ, exp_pnt)
		val = res**2
		self.assertEqual(np.unique(exp_pnt.gdf.value)[0], val)
		lat = np.array([574891.12225222, 547411.67251407, 499789.43052324,
		460177.36473906, 364182.25061015, 369862.57549558,
		394014.84676059, 444749.18166595, 514229.75466288,
		568740.1294081 , 507005.3662343 , 458457.48789088])
		np.testing.assert_allclose(
		exp_pnt.gdf.groupby(level=[0])['latitude'].nth(0).values,
		lat, atol=100000)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I really have trouble understanding this test 🙈 Can you walk me through what you are trying to test exactly?

I think one important check would be if _interp_one_poly is executed instead of _interp_one_poly_m when passing an Exposures object in a projected CRS. That one could be solved with mocking.

The final assert_allclose with an absolute tolernace of 100000 seems pointless to me. Can you explain why the tolerance is that large?

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're right, strictly speaking one should test that _interp_one_poly_m does not get called and _interp_one_poly is called instead in this example. But that surpassed my leisurely testing skills :D

I re-did this test now, just checking that the metre-based projection wasnt changed by the to_metres command (which indirectly checks this).

Not sure how mocking works, but could be an option for (future?) re-writes of the lines_polys_tests? It still has untested functions, like the ones you mentioned..

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Anyways, this is not what's happening apparently (jenkins complains that the crs was overwritten to epsg:4326) - a behaviour which i cannot reproduce locally, where all the tests pass. will try to find out how this is possible

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can come up with a code that mocks of the _interp_one_poly[_m] methods. We can also discuss that in person.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

jenkins complains that the crs was overwritten to epsg:4326

@Evelyn-M This is probably because the code uses pd.concat to concatenate an empty (EPSG:4326) GeoDataFrame with the created one, and places it into a new GeoDataFrame. Either the empty frame overwrites the CRS of the created one, or the CRS is dropped along the way by pd.concat. See

gdf_pnt = gpd.GeoDataFrame([])

and following lines

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The bug was not in that line, but in the function swapping the geometry column. It happens that when one renames the geometry column, the crs information is lost.

Fixed in 87f2e8e

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Additional info: line 525, which initializes an empty geodataframe does not fix a crs. When an empty geodataframe is created, not geometry column is defined, and thus no crs.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

great, nice work. thanks both @chahank , @peanutfun !
I guess we can merge and close this?

Co-authored-by: Lukas Riedel <34276446+peanutfun@users.noreply.github.com>

Just one small note here: the reason why this check was previously not implemented is that when you set to_meters=True, each polygon is mapped onto an equal-area projection centered on the center of the polygon. Thus, if a user has polygons all over the world, the local projections help in minimizing the area errors.

Now, a projected crs does not necessarily preserve area, and if global, leads to local errors. With the proposed change, an imho important use case is thus now not possible anymore (i.e. reprojecting polygons from a projected crs to a local projected crs with area preservation).

Then you'd need to code up something for this usecase ;) at the moment, the search for the best-fitting metre-based projection expects a lat/lon based crs to start with. Hence the fail.

Then you'd need to code up something for this usecase ;) at the moment, the search for the best-fitting metre-based projection expects a lat/lon based crs to start with. Hence the fail.

aaa yeah, right, good point.

…python into hotfix/lines_polys_metres

Evelyn-M deleted the hotfix/lines_polys_metres branch

March 3, 2023 16:37