Commit 1add7f60 authored by Clément Haëck's avatar Clément Haëck
Browse files

Add 1km to make_land

Add binary closing to land_large (with large kernel
in convolution, it makes holes a few pixels large)
parent 64f0ddf6
...@@ -3,10 +3,12 @@ ...@@ -3,10 +3,12 @@
import numpy as np import numpy as np
import xarray as xr import xarray as xr
from global_land_mask import globe
from scipy import ndimage as nd from scipy import ndimage as nd
import lib import lib
import lib.data.globcolour import lib.data.globcolour
import lib.data.modis
import lib.data.ostia import lib.data.ostia
import lib.zones import lib.zones
...@@ -15,6 +17,7 @@ def main(): ...@@ -15,6 +17,7 @@ def main():
args = lib.get_args(['region']) args = lib.get_args(['region'])
make_land(args) make_land(args)
def get_filename(grid, args): def get_filename(grid, args):
filename = lib.zones.get_filename('land_mask', args, grid=grid) filename = lib.zones.get_filename('land_mask', args, grid=grid)
print(filename) print(filename)
...@@ -25,6 +28,7 @@ def write_land(ds, land, filename, grid): ...@@ -25,6 +28,7 @@ def write_land(ds, land, filename, grid):
ds = ds.drop_vars(ds.data_vars) ds = ds.drop_vars(ds.data_vars)
ds.attrs = dict(grid=grid) ds.attrs = dict(grid=grid)
ds['land'] = land ds['land'] = land
if 'time' in ds.land.dims:
ds = ds.isel(time=0) ds = ds.isel(time=0)
ds = enlarge(ds) ds = enlarge(ds)
ds.to_netcdf(filename, encoding={'land': {'zlib': True}}) ds.to_netcdf(filename, encoding={'land': {'zlib': True}})
...@@ -45,10 +49,28 @@ def enlarge(ds): ...@@ -45,10 +49,28 @@ def enlarge(ds):
input_core_dims=[['lat', 'lon']], input_core_dims=[['lat', 'lon']],
output_core_dims=[['lat', 'lon']], output_core_dims=[['lat', 'lon']],
kwargs=dict(weights=kernel, mode='nearest')) kwargs=dict(weights=kernel, mode='nearest'))
ds['closed'] = xr.apply_ufunc(
nd.binary_closing, ds['land_large'],
input_core_dims=[['lat', 'lon']],
output_core_dims=[['lat', 'lon']]
)
ds['land_large'] = ds['land_large'] | ds['closed']
ds.drop_vars('closed')
return ds return ds
def make_land(args): def make_land(args):
# 1km
grid = '1km'
filename = get_filename(grid, args)
finder = lib.data.modis.get_finder(args)
ds = xr.open_dataset(finder.get_files()[0])
lon_grid, lat_grid = np.meshgrid(ds.lon.values, ds.lat.values)
land = globe.is_land(lat_grid, lon_grid)
ds['land'] = (['lat', 'lon'], land)
write_land(ds, ds['land'], filename, grid)
# EPSG32662 # EPSG32662
grid = '4km_EPSG32662' grid = '4km_EPSG32662'
filename = get_filename(grid, args) filename = get_filename(grid, args)
......
...@@ -28,4 +28,5 @@ root_data: /... ...@@ -28,4 +28,5 @@ root_data: /...
shapely shapely
xarray xarray
filefinder filefinder
global-land-mask
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment