-
Notifications
You must be signed in to change notification settings - Fork 127
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
Fix field plotting #1247
Fix field plotting #1247
Changes from 12 commits
0450573
7ff9c89
2aa4cc4
a46c891
33aebc4
75eac5d
bccb214
b6c924d
5cf5c2a
88e999e
28f86c0
4d2f340
d06d1a9
9ff370b
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -28,4 +28,5 @@ dependencies: | |
- nbval | ||
- scikit-learn | ||
- pykdtree | ||
- cartopy | ||
- zarr |
Original file line number | Diff line number | Diff line change | ||
---|---|---|---|---|
|
@@ -19,10 +19,11 @@ dependencies: | |||
- scipy>=0.16.0 | ||||
- tqdm | ||||
- xarray>=0.5.1 | ||||
- dask>=2.0 | ||||
- dask<=2022.9.0 | ||||
- cftime>=1.3.1 | ||||
- ipykernel | ||||
- pytest | ||||
- nbval | ||||
- pykdtree | ||||
- cartopy | ||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. See above, I suggest not to make cartopy an explicit dependency (as we had not done before)
Suggested change
|
||||
- zarr |
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -2,6 +2,7 @@ | |
from datetime import timedelta as delta | ||
|
||
import numpy as np | ||
import dask.array as da | ||
import copy | ||
|
||
from parcels.field import Field | ||
|
@@ -101,7 +102,7 @@ def plotparticles(particles, with_particles=True, show_time=None, field=None, do | |
|
||
|
||
def plotfield(field, show_time=None, domain=None, depth_level=0, projection='PlateCarree', land=True, | ||
vmin=None, vmax=None, savefile=None, **kwargs): | ||
vmin=None, vmax=None, savefile=None, use3D=False, **kwargs): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Is this variable still needed now? I am a bit confused what this does There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. possibly not, indeed |
||
"""Function to plot a Parcels Field | ||
|
||
:param show_time: Time in seconds from start after which to show the Field | ||
|
@@ -113,6 +114,7 @@ def plotfield(field, show_time=None, domain=None, depth_level=0, projection='Pla | |
:param vmax: maximum colour scale (only in single-plot mode) | ||
:param savefile: Name of a file to save the plot to | ||
:param animation: Boolean whether result is a single plot, or an animation | ||
:param use3D: tells if requested projection is 2D or 3D | ||
""" | ||
|
||
if type(field) is VectorField: | ||
|
@@ -130,7 +132,7 @@ def plotfield(field, show_time=None, domain=None, depth_level=0, projection='Pla | |
logger.warning('Field.show() does not always correctly determine the domain for curvilinear grids. ' | ||
'Use plotting with caution and perhaps use domain argument as in the NEMO 3D tutorial') | ||
|
||
plt, fig, ax, cartopy = create_parcelsfig_axis(spherical, land, projection=projection, cartopy_features=kwargs.pop('cartopy_features', [])) | ||
plt, fig, ax, cartopy = create_parcelsfig_axis(spherical, land, projection=projection, cartopy_features=kwargs.pop('cartopy_features', []), use3D=use3D) | ||
if plt is None: | ||
return None, None, None, None # creating axes was not possible | ||
|
||
|
@@ -165,6 +167,8 @@ def plotfield(field, show_time=None, domain=None, depth_level=0, projection='Pla | |
data[i] = np.squeeze(fld.data)[depth_level, latS:latN, lonW:lonE] | ||
else: | ||
data[i] = np.squeeze(fld.data)[latS:latN, lonW:lonE] | ||
if isinstance(data[i], da.Array): | ||
data[i] = np.array(data[i]) | ||
|
||
if plottype == 'vector': | ||
if field[0].interp_method == 'cgrid_velocity': | ||
|
@@ -264,21 +268,26 @@ def plotfield(field, show_time=None, domain=None, depth_level=0, projection='Pla | |
return plt, fig, ax, cartopy | ||
|
||
|
||
def create_parcelsfig_axis(spherical, land=True, projection='PlateCarree', central_longitude=0, cartopy_features=[]): | ||
def create_parcelsfig_axis(spherical, land=True, projection='PlateCarree', central_longitude=0, cartopy_features=[], use3D=False): | ||
try: | ||
import matplotlib.pyplot as plt | ||
except: | ||
logger.info("Visualisation is not possible. Matplotlib not found.") | ||
return None, None, None, None # creating axes was not possible | ||
|
||
if spherical and projection: | ||
if use3D: | ||
cartopy = None | ||
fig, ax = plt.subplots(1, 1, subplot_kw={'projection': '3d'}) | ||
ax.grid() | ||
elif spherical and projection: | ||
try: | ||
import cartopy | ||
except: | ||
logger.info("Visualisation of field with geographic coordinates is not possible. Cartopy not found.") | ||
return None, None, None, None # creating axes was not possible | ||
|
||
projection = cartopy.crs.PlateCarree(central_longitude) if projection == 'PlateCarree' else projection | ||
# projection = '3d' if use3D else projection | ||
CKehl marked this conversation as resolved.
Show resolved
Hide resolved
|
||
fig, ax = plt.subplots(1, 1, subplot_kw={'projection': projection}) | ||
try: # gridlines not supported for all projections | ||
if isinstance(projection, cartopy.crs.PlateCarree) and central_longitude != 0: | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,54 @@ | ||
from parcels import FieldSet, ParticleSet, JITParticle, AdvectionRK4, ErrorCode | ||
from datetime import timedelta as delta | ||
import numpy as np | ||
import dask.array as da # NOQA | ||
from os import path | ||
from matplotlib import pyplot as plt # NOQA | ||
import pytest | ||
|
||
|
||
def periodicBC(particle, fieldset, time): | ||
if particle.lon > 180: | ||
particle.lon -= 360 | ||
|
||
|
||
def test_field_from_netcdf(): | ||
data_path = path.join(path.dirname(__file__), 'test_data/') | ||
|
||
filenames = {'U': {'lon': data_path + 'mask_nemo_cross_180lon.nc', | ||
'lat': data_path + 'mask_nemo_cross_180lon.nc', | ||
'data': data_path + 'Uu_eastward_nemo_cross_180lon.nc'}, | ||
'V': {'lon': data_path + 'mask_nemo_cross_180lon.nc', | ||
'lat': data_path + 'mask_nemo_cross_180lon.nc', | ||
'data': data_path + 'Vv_eastward_nemo_cross_180lon.nc'} | ||
} | ||
variables = {'U': 'U', | ||
'V': 'V'} | ||
dimensions = {'lon': 'glamf', 'lat': 'gphif'} | ||
return FieldSet.from_netcdf(filenames, variables, dimensions, interp_method='cgrid_velocity', chunksize='auto', allow_time_extrapolation=True) | ||
|
||
|
||
@pytest.fixture(name="fieldset") | ||
def fieldset_fixture(): | ||
return test_field_from_netcdf() | ||
|
||
|
||
def test_pset_create_field(fieldset, npart=100): | ||
lonp = -180 * np.ones(npart) | ||
latp = [i for i in np.linspace(-70, 88, npart)] | ||
pset = ParticleSet.from_list(fieldset, JITParticle, lon=lonp, lat=latp) | ||
return pset | ||
|
||
|
||
def DeleteParticle(particle, fieldset, time): | ||
particle.delete() | ||
|
||
|
||
if __name__ == '__main__': | ||
fset = test_field_from_netcdf() | ||
print(fset) | ||
pset = test_pset_create_field(fset) | ||
kernels = pset.Kernel(AdvectionRK4) + periodicBC | ||
pset.execute(kernels, dt=delta(hours=1), output_file=None, | ||
recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle}) | ||
pset.show(field=fset.U) | ||
Comment on lines
+1
to
+54
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. What does this new set of functions test that was not tested before? Why is it needed? |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why is cartopy now need? That was not the case before, and I suggest not to keep it
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
that will make the new test fail that actually tests platecarree plotting