Skip to content

Commit

Permalink
Use equal area Mollweide projection for computing fraction of overlap…
Browse files Browse the repository at this point in the history
… between inventory and model grid cells (#4)
  • Loading branch information
jmhaussaire authored Apr 1, 2020
1 parent 1f3528c commit 7628b63
Showing 1 changed file with 6 additions and 5 deletions.
11 changes: 6 additions & 5 deletions emiproc/grids.py
Original file line number Diff line number Diff line change
Expand Up @@ -743,18 +743,19 @@ def intersected_cells(self, corners):
# The inventory cell lies outside the cosmo grid
return []

# Here we assume a flat earth. The error is less than 1% for typical
# grid sizes over europe. Since we're interested in the ratio of areas,
# we can calculate in degrees^2.
molly = ccrs.Mollweide()
corners = molly.transform_points(self.projection,corners[:,0],corners[:,1])
inv_cell = Polygon(corners)


intersections = []
# make sure we iterate only over valid gridpoint indices
for i in range(max(0, lon_idx_min), min(self.nx, lon_idx_max)):
for j in range(max(0, lat_idx_min), min(self.ny, lat_idx_max)):
corners = list(zip(*self.cell_corners(i, j)))
corners = np.array(list(zip(*self.cell_corners(i, j))))
corners = molly.transform_points(self.projection,corners[:,0],corners[:,1])

cosmo_cell = Polygon(corners)

if cosmo_cell.intersects(inv_cell):
overlap = cosmo_cell.intersection(inv_cell)
intersections.append((i, j, overlap.area / inv_cell.area))
Expand Down

0 comments on commit 7628b63

Please sign in to comment.