Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Don't re-project projected polygon gdfs in metre-disaggregation #666

Merged
merged 11 commits into from
Mar 3, 2023

Conversation

Evelyn-M
Copy link
Collaborator

@Evelyn-M Evelyn-M commented Feb 28, 2023

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

PR Reviewer Checklist

Copy link
Member

@peanutfun peanutfun left a comment

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

climada/util/lines_polys_handler.py Outdated Show resolved Hide resolved
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)
Copy link
Member

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?

Copy link
Collaborator Author

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..

Copy link
Collaborator Author

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

Copy link
Member

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.

Copy link
Member

Choose a reason for hiding this comment

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

@Evelyn-M done!

Copy link
Member

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

Copy link
Member

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

Copy link
Member

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.

Copy link
Collaborator Author

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?

@chahank
Copy link
Member

chahank commented Mar 1, 2023

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).

@Evelyn-M
Copy link
Collaborator Author

Evelyn-M commented Mar 1, 2023

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.

@chahank
Copy link
Member

chahank commented Mar 1, 2023

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants