From 20c3ba3a633ec165e200a5e2201e31d792d16583 Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Thu, 1 Apr 2021 16:45:59 +0100 Subject: [PATCH 1/9] add new module for plotting UM domains #188 --- lib/catnip/um_domain.py | 204 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 204 insertions(+) create mode 100644 lib/catnip/um_domain.py diff --git a/lib/catnip/um_domain.py b/lib/catnip/um_domain.py new file mode 100644 index 0000000..48e544e --- /dev/null +++ b/lib/catnip/um_domain.py @@ -0,0 +1,204 @@ +''' +Plot the LAM domains +Variable resolution not supported. Old dynamics +is supported. +''' +import cartopy +import cartopy.crs as ccrs +import cartopy.feature as cfeature +import iris +import matplotlib.pyplot as plt +import numpy as np +import os + +class UmDomain: + ''' + A class to hold the information needed + to define a model domain in the MetUM + ''' + + def __init__(self, delta_lat, delta_lon, + frstlata, frstlona, + polelata, polelona, + global_row_length, global_rows, + old_dynam=False): + + 'Initialise the class holding a rotated pole domain' + self.delta_lat = delta_lat # lattitude grid spacing + self.delta_lon = delta_lon # longitude grid spacing + self.frstlata = frstlata # first latitude + self.frstlona = frstlona # first lon + self.polelata = polelata # latidude of rotated pole + self.polelona = polelona # longitude of rotated pole + self.global_row_length = global_row_length # no of columns + self.global_rows = global_rows # no of rows + self.old_dynam = old_dynam # is this an old dynamics grid? + + self.theta_cube_2d = None + self.theta_cube_3d = None + + def last_lon_lat(self): + 'get the last lat and lon in the grid' + # NB is the right for old dynamics? Only matters + # if going to use it to make old dynamics cubes + + lastlona = (self.frstlona + + (self.global_row_length - 1) * self.delta_lon) + lastlata = (self.frstlata + + (self.global_rows - 1) * self.delta_lat) + + return lastlona, lastlata + + def make_2d_theta_cube(self): + ''' + using the information on the grid, make a 2D cube + of theta + ''' + + # calculate the lat and lon arrays + lastlona, lastlata = self.last_lon_lat() + theta_lon = np.linspace(self.frstlona, lastlona, + num=self.global_row_length) + theta_lat = np.linspace(self.frstlata, lastlata, + num=self.global_rows) + + # create a 2D numpy ndarray + theta_data = np.zeros((len(theta_lon), len(theta_lat))) + + lon_coord = iris.coords.DimCoord(theta_lon, + standard_name='grid_longitude', + units='degrees', + coord_system= + iris.coord_systems.RotatedGeogCS(self.polelata, self.polelona)) + lat_coord = iris.coords.DimCoord(theta_lat, + standard_name='grid_latitude', + units='degrees', + coord_system= + iris.coord_systems.RotatedGeogCS(self.polelata, self.polelona) + ) + + self.theta_cube_2d = iris.cube.Cube(theta_data, long_name = 'air_potential_temperature', + units = 'K', + dim_coords_and_dims=[ + (lon_coord, 0), + (lat_coord, 1)]) + + return self.theta_cube_2d + + def get_corners(self): + ''' + Find the lat/lon of the extreme corners of the + domain. + Note that for old dynamics first lat and lon + are the grid box centres whereas for ND/EG + they are the grid box corners + ''' + + # Old dynamics + if self.old_dynam: + + # add code here + # longitude - as this refers to centres, + # we have to subtract half a grid box + self.minlon = self.frstlona - self.delta_lon/2. + # latitude - not that frstlona is + # the maximum latitude + self.minlat = (self.frstlata + - self.delta_lon * (self.global_rows - 0.5)) + + # max lon/lat + self.maxlon = (self.frstlona + + (self.global_row_length -0.5) * self.delta_lon) + self.maxlat = self.frstlata + self.delta_lat/2. + + + # New dynamics/ENDGAME + else: + + # min lon and min lat are same as firstlona etc + self.minlon = self.frstlona + self.minlat = self.frstlata + + # as this refers to grid box corners not + # centres we multiply by N not by N-1 + self.maxlat = (self.frstlata + + self.global_rows * self.delta_lat) + + self.maxlon = (self.frstlona + + self.global_row_length * self.delta_lon) + + + return self.minlat, self.minlon, self.maxlat, self.maxlon + +def domain_plotter(title, mapextent, domains, background=None, resolution='1.0', + central_longitude=0.,plot_origin=True): + ''' + Take the standard information used to describe + domains for the UM and plot them as a boxes on a map + + Could add default extent as whole globe? + + extent - left, right, bottom, top in plate carree coordinates + ''' + + # create axes, with lower threshold for interpolation + # so that the lines are smoothly curving + pc = ccrs.PlateCarree(central_longitude=central_longitude) + pc._threshold /= 10 + + ax = plt.axes(projection=pc) + ax.set_extent(mapextent) + if background is None: + ax.stock_img() + else: + # tell cartopy where to find backgrounds + os.environ["CARTOPY_USER_BACKGROUNDS"] = '/data/users/fris/nasa_bg' + ax.background_img(name=background, resolution=resolution) + ax.coastlines(color='gray') + # Create a feature for States/Admin 1 regions at 1:50m from Natural Earth + states_provinces = cfeature.NaturalEarthFeature( + category='cultural', + name='admin_1_states_provinces_lines', + scale='10m', + facecolor='none') + + ax.add_feature(states_provinces, edgecolor='gray') + ax.add_feature(cartopy.feature.RIVERS) + ax.add_feature(cartopy.feature.BORDERS) + ax.add_feature(cartopy.feature.LAKES) + ax.add_feature(cartopy.feature.COASTLINE) + gl = ax.gridlines(draw_labels=True) + gl.xlabels_top = False + + + plt.title(title) + + for domain in domains: + rotated_pole = ccrs.RotatedPole(pole_latitude=domain.polelata, + pole_longitude=domain.polelona) + + # specify x and y of corners in rotated pole co-ordinates, + # returning to start to make a closed curve + + # Replace this with a call to domain.get_corners() + # and then + # x = domain.minlon, domain.minlon, domain.maxlon, + # domain.maxlon, domain.minlon] + + domain.get_corners() + + x = [domain.minlon, domain.minlon, domain.maxlon, + domain.maxlon, domain.minlon] + y = [domain.minlat, domain.maxlat, domain.maxlat, + domain.minlat, domain.minlat] + + + p = ax.plot(x, y, marker='o', transform=rotated_pole) + color = p[0].get_color() + if plot_origin: + # Add a marker for the origin of the domain + ax.plot(0.0, 0.0, marker='o', color=color, markersize=5, + alpha=0.7, transform=rotated_pole) + + plt.show() + plt.close() From 1f40646f32940a29e9e1adb79e3836e53a96e9a1 Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Thu, 8 Apr 2021 14:59:02 +0100 Subject: [PATCH 2/9] move calls to plot to separate routine --- lib/catnip/test_plotter.py | 73 ++++++++++++++++++++++++++++++++++++++ lib/catnip/um_domain.py | 3 +- 2 files changed, 74 insertions(+), 2 deletions(-) create mode 100644 lib/catnip/test_plotter.py diff --git a/lib/catnip/test_plotter.py b/lib/catnip/test_plotter.py new file mode 100644 index 0000000..a70178d --- /dev/null +++ b/lib/catnip/test_plotter.py @@ -0,0 +1,73 @@ +# Test plotter code +# may not go into final commit but useful here for now +import matplotlib.pyplot as plt +import cartopy.io.shapereader as shpreader +import cartopy.crs as ccrs +import um_domain + +pole_lat = 70.0 +pole_lon = 310.0 + +# 12lm domain-24 +domain_12k = um_domain.UmDomain(delta_lat=0.11, delta_lon=0.11, + frstlata=-24.0, + frstlona=328.0, + polelata=pole_lat, polelona=pole_lon, + global_row_length=700, global_rows=550) + +domain_4k_v2 = um_domain.UmDomain(delta_lat=0.0405, delta_lon=0.0405, + frstlata=0., frstlona=343.0, + polelata=pole_lat, polelona=pole_lon, + global_row_length=500, global_rows=432) + +# plot all of these on one map + +mapextent = [80, 195, -10, 45] # x1, x2, y1, y2 + +title = '4 & 12 km domains. ' +fig = um_domain.domain_plotter(title, mapextent, [domain_4k_v2, domain_12k], + background='topo', resolution='0.25', + central_longitude=180) + +# list of places to add to the plot +site_list = ['Beijing', 'Shanghai', 'Guangzhou', 'Shenzhen', 'Xian', + 'Chongqing', 'Chengdu', 'Harbin', 'Wuhan'] + +# load populated places data from natural Earth +shpfilename = shpreader.natural_earth(resolution='10m', category='cultural', + name='populated_places') +reader = shpreader.Reader(shpfilename) +places = reader.records() + +# filter only sites in the list above +sites_shapes = [place for place in places if place.attributes['NAME'] in site_list ] + +ax = plt.gca() + +ax.scatter([point.attributes['LONGITUDE'] for point in sites_shapes], + [point.attributes['LATITUDE'] for point in sites_shapes], + c='red', + transform=ccrs.Geodetic()) +plt.show() +plt.close() + +# and again for 4 km and inner sites +site_list = ['Shanghai', 'Guangzhou', 'Shenzhen', 'Wuhan'] + +# filter only sites in the list above +places = reader.records() +sites_shapes = [place for place in places if place.attributes['NAME'] in site_list ] +print(len(sites_shapes)) + +mapextent = [108, 136, 18, 40] # x1, x2, y1, y2 +title = '4 km domain' +fig = um_domain.domain_plotter(title, mapextent, [domain_4k_v2], + background='topo', resolution='0.25', + central_longitude=180) +ax = plt.gca() + +ax.scatter([point.attributes['LONGITUDE'] for point in sites_shapes], + [point.attributes['LATITUDE'] for point in sites_shapes], + c='red', transform=ccrs.Geodetic()) +plt.show() +plt.close() diff --git a/lib/catnip/um_domain.py b/lib/catnip/um_domain.py index 48e544e..75a5e00 100644 --- a/lib/catnip/um_domain.py +++ b/lib/catnip/um_domain.py @@ -200,5 +200,4 @@ def domain_plotter(title, mapextent, domains, background=None, resolution='1.0', ax.plot(0.0, 0.0, marker='o', color=color, markersize=5, alpha=0.7, transform=rotated_pole) - plt.show() - plt.close() + return plt.gcf() From 06c00c9b0efe041107d849e784c60ec0dafbb847 Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Thu, 8 Apr 2021 15:40:34 +0100 Subject: [PATCH 3/9] move code to plot cities on a map to own function --- lib/catnip/test_plotter.py | 26 ++------------------------ lib/catnip/um_domain.py | 23 +++++++++++++++++++++++ 2 files changed, 25 insertions(+), 24 deletions(-) diff --git a/lib/catnip/test_plotter.py b/lib/catnip/test_plotter.py index a70178d..5807dc6 100644 --- a/lib/catnip/test_plotter.py +++ b/lib/catnip/test_plotter.py @@ -33,41 +33,19 @@ site_list = ['Beijing', 'Shanghai', 'Guangzhou', 'Shenzhen', 'Xian', 'Chongqing', 'Chengdu', 'Harbin', 'Wuhan'] -# load populated places data from natural Earth -shpfilename = shpreader.natural_earth(resolution='10m', category='cultural', - name='populated_places') -reader = shpreader.Reader(shpfilename) -places = reader.records() - -# filter only sites in the list above -sites_shapes = [place for place in places if place.attributes['NAME'] in site_list ] - -ax = plt.gca() - -ax.scatter([point.attributes['LONGITUDE'] for point in sites_shapes], - [point.attributes['LATITUDE'] for point in sites_shapes], - c='red', - transform=ccrs.Geodetic()) +um_domain.plot_cities(site_list) plt.show() plt.close() # and again for 4 km and inner sites site_list = ['Shanghai', 'Guangzhou', 'Shenzhen', 'Wuhan'] -# filter only sites in the list above -places = reader.records() -sites_shapes = [place for place in places if place.attributes['NAME'] in site_list ] -print(len(sites_shapes)) - mapextent = [108, 136, 18, 40] # x1, x2, y1, y2 title = '4 km domain' fig = um_domain.domain_plotter(title, mapextent, [domain_4k_v2], background='topo', resolution='0.25', central_longitude=180) -ax = plt.gca() +um_domain.plot_cities(site_list) -ax.scatter([point.attributes['LONGITUDE'] for point in sites_shapes], - [point.attributes['LATITUDE'] for point in sites_shapes], - c='red', transform=ccrs.Geodetic()) plt.show() plt.close() diff --git a/lib/catnip/um_domain.py b/lib/catnip/um_domain.py index 75a5e00..66115ca 100644 --- a/lib/catnip/um_domain.py +++ b/lib/catnip/um_domain.py @@ -6,6 +6,7 @@ import cartopy import cartopy.crs as ccrs import cartopy.feature as cfeature +import cartopy.io.shapereader as shpreader import iris import matplotlib.pyplot as plt import numpy as np @@ -201,3 +202,25 @@ def domain_plotter(title, mapextent, domains, background=None, resolution='1.0', alpha=0.7, transform=rotated_pole) return plt.gcf() + +def plot_cities(site_list): + ''' + Given a list of names of cities, plot them on the current axes + Must be included in the Natural Earth populated places + ''' + + # load populated places data from natural Earth + shpfilename = shpreader.natural_earth(resolution='10m', category='cultural', + name='populated_places') + reader = shpreader.Reader(shpfilename) + places = reader.records() + + # filter only sites in the list above + sites_shapes = [place for place in places if place.attributes['NAME'] in site_list ] + + ax = plt.gca() + + ax.scatter([point.attributes['LONGITUDE'] for point in sites_shapes], + [point.attributes['LATITUDE'] for point in sites_shapes], + c='red', + transform=ccrs.Geodetic()) From 783f8b8da9be6f93dc4d1d275fd557cda500b110 Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Mon, 16 Aug 2021 16:08:07 +0100 Subject: [PATCH 4/9] modify to be the new proposals for CNCP runs --- lib/catnip/test_plotter.py | 12 ++++++------ lib/catnip/um_domain.py | 4 ++-- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/lib/catnip/test_plotter.py b/lib/catnip/test_plotter.py index 5807dc6..5115f73 100644 --- a/lib/catnip/test_plotter.py +++ b/lib/catnip/test_plotter.py @@ -11,20 +11,20 @@ # 12lm domain-24 domain_12k = um_domain.UmDomain(delta_lat=0.11, delta_lon=0.11, frstlata=-24.0, - frstlona=328.0, + frstlona=310.0, polelata=pole_lat, polelona=pole_lon, global_row_length=700, global_rows=550) domain_4k_v2 = um_domain.UmDomain(delta_lat=0.0405, delta_lon=0.0405, - frstlata=0., frstlona=343.0, + frstlata=-2., frstlona=340.0, polelata=pole_lat, polelona=pole_lon, global_row_length=500, global_rows=432) # plot all of these on one map -mapextent = [80, 195, -10, 45] # x1, x2, y1, y2 +mapextent = [50, 178, -10, 45] # x1, x2, y1, y2 -title = '4 & 12 km domains. ' +title = ' New 4 & 12 km domains. ' fig = um_domain.domain_plotter(title, mapextent, [domain_4k_v2, domain_12k], background='topo', resolution='0.25', central_longitude=180) @@ -40,8 +40,8 @@ # and again for 4 km and inner sites site_list = ['Shanghai', 'Guangzhou', 'Shenzhen', 'Wuhan'] -mapextent = [108, 136, 18, 40] # x1, x2, y1, y2 -title = '4 km domain' +mapextent = [104, 132, 12, 36] # x1, x2, y1, y2 +title = 'New 4 km domain' fig = um_domain.domain_plotter(title, mapextent, [domain_4k_v2], background='topo', resolution='0.25', central_longitude=180) diff --git a/lib/catnip/um_domain.py b/lib/catnip/um_domain.py index 66115ca..4a35ff7 100644 --- a/lib/catnip/um_domain.py +++ b/lib/catnip/um_domain.py @@ -194,7 +194,7 @@ def domain_plotter(title, mapextent, domains, background=None, resolution='1.0', domain.minlat, domain.minlat] - p = ax.plot(x, y, marker='o', transform=rotated_pole) + p = ax.plot(x, y, marker=None, transform=rotated_pole) color = p[0].get_color() if plot_origin: # Add a marker for the origin of the domain @@ -206,7 +206,7 @@ def domain_plotter(title, mapextent, domains, background=None, resolution='1.0', def plot_cities(site_list): ''' Given a list of names of cities, plot them on the current axes - Must be included in the Natural Earth populated places + Must be included in the Natural Earth populated places ''' # load populated places data from natural Earth From f5dfbe427821dad64291f3b98d740916fe6d57fa Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Thu, 1 Apr 2021 16:45:59 +0100 Subject: [PATCH 5/9] add new module for plotting UM domains #188 --- lib/catnip/um_domain.py | 204 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 204 insertions(+) create mode 100644 lib/catnip/um_domain.py diff --git a/lib/catnip/um_domain.py b/lib/catnip/um_domain.py new file mode 100644 index 0000000..48e544e --- /dev/null +++ b/lib/catnip/um_domain.py @@ -0,0 +1,204 @@ +''' +Plot the LAM domains +Variable resolution not supported. Old dynamics +is supported. +''' +import cartopy +import cartopy.crs as ccrs +import cartopy.feature as cfeature +import iris +import matplotlib.pyplot as plt +import numpy as np +import os + +class UmDomain: + ''' + A class to hold the information needed + to define a model domain in the MetUM + ''' + + def __init__(self, delta_lat, delta_lon, + frstlata, frstlona, + polelata, polelona, + global_row_length, global_rows, + old_dynam=False): + + 'Initialise the class holding a rotated pole domain' + self.delta_lat = delta_lat # lattitude grid spacing + self.delta_lon = delta_lon # longitude grid spacing + self.frstlata = frstlata # first latitude + self.frstlona = frstlona # first lon + self.polelata = polelata # latidude of rotated pole + self.polelona = polelona # longitude of rotated pole + self.global_row_length = global_row_length # no of columns + self.global_rows = global_rows # no of rows + self.old_dynam = old_dynam # is this an old dynamics grid? + + self.theta_cube_2d = None + self.theta_cube_3d = None + + def last_lon_lat(self): + 'get the last lat and lon in the grid' + # NB is the right for old dynamics? Only matters + # if going to use it to make old dynamics cubes + + lastlona = (self.frstlona + + (self.global_row_length - 1) * self.delta_lon) + lastlata = (self.frstlata + + (self.global_rows - 1) * self.delta_lat) + + return lastlona, lastlata + + def make_2d_theta_cube(self): + ''' + using the information on the grid, make a 2D cube + of theta + ''' + + # calculate the lat and lon arrays + lastlona, lastlata = self.last_lon_lat() + theta_lon = np.linspace(self.frstlona, lastlona, + num=self.global_row_length) + theta_lat = np.linspace(self.frstlata, lastlata, + num=self.global_rows) + + # create a 2D numpy ndarray + theta_data = np.zeros((len(theta_lon), len(theta_lat))) + + lon_coord = iris.coords.DimCoord(theta_lon, + standard_name='grid_longitude', + units='degrees', + coord_system= + iris.coord_systems.RotatedGeogCS(self.polelata, self.polelona)) + lat_coord = iris.coords.DimCoord(theta_lat, + standard_name='grid_latitude', + units='degrees', + coord_system= + iris.coord_systems.RotatedGeogCS(self.polelata, self.polelona) + ) + + self.theta_cube_2d = iris.cube.Cube(theta_data, long_name = 'air_potential_temperature', + units = 'K', + dim_coords_and_dims=[ + (lon_coord, 0), + (lat_coord, 1)]) + + return self.theta_cube_2d + + def get_corners(self): + ''' + Find the lat/lon of the extreme corners of the + domain. + Note that for old dynamics first lat and lon + are the grid box centres whereas for ND/EG + they are the grid box corners + ''' + + # Old dynamics + if self.old_dynam: + + # add code here + # longitude - as this refers to centres, + # we have to subtract half a grid box + self.minlon = self.frstlona - self.delta_lon/2. + # latitude - not that frstlona is + # the maximum latitude + self.minlat = (self.frstlata + - self.delta_lon * (self.global_rows - 0.5)) + + # max lon/lat + self.maxlon = (self.frstlona + + (self.global_row_length -0.5) * self.delta_lon) + self.maxlat = self.frstlata + self.delta_lat/2. + + + # New dynamics/ENDGAME + else: + + # min lon and min lat are same as firstlona etc + self.minlon = self.frstlona + self.minlat = self.frstlata + + # as this refers to grid box corners not + # centres we multiply by N not by N-1 + self.maxlat = (self.frstlata + + self.global_rows * self.delta_lat) + + self.maxlon = (self.frstlona + + self.global_row_length * self.delta_lon) + + + return self.minlat, self.minlon, self.maxlat, self.maxlon + +def domain_plotter(title, mapextent, domains, background=None, resolution='1.0', + central_longitude=0.,plot_origin=True): + ''' + Take the standard information used to describe + domains for the UM and plot them as a boxes on a map + + Could add default extent as whole globe? + + extent - left, right, bottom, top in plate carree coordinates + ''' + + # create axes, with lower threshold for interpolation + # so that the lines are smoothly curving + pc = ccrs.PlateCarree(central_longitude=central_longitude) + pc._threshold /= 10 + + ax = plt.axes(projection=pc) + ax.set_extent(mapextent) + if background is None: + ax.stock_img() + else: + # tell cartopy where to find backgrounds + os.environ["CARTOPY_USER_BACKGROUNDS"] = '/data/users/fris/nasa_bg' + ax.background_img(name=background, resolution=resolution) + ax.coastlines(color='gray') + # Create a feature for States/Admin 1 regions at 1:50m from Natural Earth + states_provinces = cfeature.NaturalEarthFeature( + category='cultural', + name='admin_1_states_provinces_lines', + scale='10m', + facecolor='none') + + ax.add_feature(states_provinces, edgecolor='gray') + ax.add_feature(cartopy.feature.RIVERS) + ax.add_feature(cartopy.feature.BORDERS) + ax.add_feature(cartopy.feature.LAKES) + ax.add_feature(cartopy.feature.COASTLINE) + gl = ax.gridlines(draw_labels=True) + gl.xlabels_top = False + + + plt.title(title) + + for domain in domains: + rotated_pole = ccrs.RotatedPole(pole_latitude=domain.polelata, + pole_longitude=domain.polelona) + + # specify x and y of corners in rotated pole co-ordinates, + # returning to start to make a closed curve + + # Replace this with a call to domain.get_corners() + # and then + # x = domain.minlon, domain.minlon, domain.maxlon, + # domain.maxlon, domain.minlon] + + domain.get_corners() + + x = [domain.minlon, domain.minlon, domain.maxlon, + domain.maxlon, domain.minlon] + y = [domain.minlat, domain.maxlat, domain.maxlat, + domain.minlat, domain.minlat] + + + p = ax.plot(x, y, marker='o', transform=rotated_pole) + color = p[0].get_color() + if plot_origin: + # Add a marker for the origin of the domain + ax.plot(0.0, 0.0, marker='o', color=color, markersize=5, + alpha=0.7, transform=rotated_pole) + + plt.show() + plt.close() From 3c0dee06fc5c2236559129c9ee054c308c8e6766 Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Thu, 8 Apr 2021 14:59:02 +0100 Subject: [PATCH 6/9] move calls to plot to separate routine --- lib/catnip/test_plotter.py | 73 ++++++++++++++++++++++++++++++++++++++ lib/catnip/um_domain.py | 3 +- 2 files changed, 74 insertions(+), 2 deletions(-) create mode 100644 lib/catnip/test_plotter.py diff --git a/lib/catnip/test_plotter.py b/lib/catnip/test_plotter.py new file mode 100644 index 0000000..a70178d --- /dev/null +++ b/lib/catnip/test_plotter.py @@ -0,0 +1,73 @@ +# Test plotter code +# may not go into final commit but useful here for now +import matplotlib.pyplot as plt +import cartopy.io.shapereader as shpreader +import cartopy.crs as ccrs +import um_domain + +pole_lat = 70.0 +pole_lon = 310.0 + +# 12lm domain-24 +domain_12k = um_domain.UmDomain(delta_lat=0.11, delta_lon=0.11, + frstlata=-24.0, + frstlona=328.0, + polelata=pole_lat, polelona=pole_lon, + global_row_length=700, global_rows=550) + +domain_4k_v2 = um_domain.UmDomain(delta_lat=0.0405, delta_lon=0.0405, + frstlata=0., frstlona=343.0, + polelata=pole_lat, polelona=pole_lon, + global_row_length=500, global_rows=432) + +# plot all of these on one map + +mapextent = [80, 195, -10, 45] # x1, x2, y1, y2 + +title = '4 & 12 km domains. ' +fig = um_domain.domain_plotter(title, mapextent, [domain_4k_v2, domain_12k], + background='topo', resolution='0.25', + central_longitude=180) + +# list of places to add to the plot +site_list = ['Beijing', 'Shanghai', 'Guangzhou', 'Shenzhen', 'Xian', + 'Chongqing', 'Chengdu', 'Harbin', 'Wuhan'] + +# load populated places data from natural Earth +shpfilename = shpreader.natural_earth(resolution='10m', category='cultural', + name='populated_places') +reader = shpreader.Reader(shpfilename) +places = reader.records() + +# filter only sites in the list above +sites_shapes = [place for place in places if place.attributes['NAME'] in site_list ] + +ax = plt.gca() + +ax.scatter([point.attributes['LONGITUDE'] for point in sites_shapes], + [point.attributes['LATITUDE'] for point in sites_shapes], + c='red', + transform=ccrs.Geodetic()) +plt.show() +plt.close() + +# and again for 4 km and inner sites +site_list = ['Shanghai', 'Guangzhou', 'Shenzhen', 'Wuhan'] + +# filter only sites in the list above +places = reader.records() +sites_shapes = [place for place in places if place.attributes['NAME'] in site_list ] +print(len(sites_shapes)) + +mapextent = [108, 136, 18, 40] # x1, x2, y1, y2 +title = '4 km domain' +fig = um_domain.domain_plotter(title, mapextent, [domain_4k_v2], + background='topo', resolution='0.25', + central_longitude=180) +ax = plt.gca() + +ax.scatter([point.attributes['LONGITUDE'] for point in sites_shapes], + [point.attributes['LATITUDE'] for point in sites_shapes], + c='red', transform=ccrs.Geodetic()) +plt.show() +plt.close() diff --git a/lib/catnip/um_domain.py b/lib/catnip/um_domain.py index 48e544e..75a5e00 100644 --- a/lib/catnip/um_domain.py +++ b/lib/catnip/um_domain.py @@ -200,5 +200,4 @@ def domain_plotter(title, mapextent, domains, background=None, resolution='1.0', ax.plot(0.0, 0.0, marker='o', color=color, markersize=5, alpha=0.7, transform=rotated_pole) - plt.show() - plt.close() + return plt.gcf() From 53b12010b6950804761af130df3f4d76fe130787 Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Thu, 8 Apr 2021 15:40:34 +0100 Subject: [PATCH 7/9] move code to plot cities on a map to own function --- lib/catnip/test_plotter.py | 26 ++------------------------ lib/catnip/um_domain.py | 23 +++++++++++++++++++++++ 2 files changed, 25 insertions(+), 24 deletions(-) diff --git a/lib/catnip/test_plotter.py b/lib/catnip/test_plotter.py index a70178d..5807dc6 100644 --- a/lib/catnip/test_plotter.py +++ b/lib/catnip/test_plotter.py @@ -33,41 +33,19 @@ site_list = ['Beijing', 'Shanghai', 'Guangzhou', 'Shenzhen', 'Xian', 'Chongqing', 'Chengdu', 'Harbin', 'Wuhan'] -# load populated places data from natural Earth -shpfilename = shpreader.natural_earth(resolution='10m', category='cultural', - name='populated_places') -reader = shpreader.Reader(shpfilename) -places = reader.records() - -# filter only sites in the list above -sites_shapes = [place for place in places if place.attributes['NAME'] in site_list ] - -ax = plt.gca() - -ax.scatter([point.attributes['LONGITUDE'] for point in sites_shapes], - [point.attributes['LATITUDE'] for point in sites_shapes], - c='red', - transform=ccrs.Geodetic()) +um_domain.plot_cities(site_list) plt.show() plt.close() # and again for 4 km and inner sites site_list = ['Shanghai', 'Guangzhou', 'Shenzhen', 'Wuhan'] -# filter only sites in the list above -places = reader.records() -sites_shapes = [place for place in places if place.attributes['NAME'] in site_list ] -print(len(sites_shapes)) - mapextent = [108, 136, 18, 40] # x1, x2, y1, y2 title = '4 km domain' fig = um_domain.domain_plotter(title, mapextent, [domain_4k_v2], background='topo', resolution='0.25', central_longitude=180) -ax = plt.gca() +um_domain.plot_cities(site_list) -ax.scatter([point.attributes['LONGITUDE'] for point in sites_shapes], - [point.attributes['LATITUDE'] for point in sites_shapes], - c='red', transform=ccrs.Geodetic()) plt.show() plt.close() diff --git a/lib/catnip/um_domain.py b/lib/catnip/um_domain.py index 75a5e00..66115ca 100644 --- a/lib/catnip/um_domain.py +++ b/lib/catnip/um_domain.py @@ -6,6 +6,7 @@ import cartopy import cartopy.crs as ccrs import cartopy.feature as cfeature +import cartopy.io.shapereader as shpreader import iris import matplotlib.pyplot as plt import numpy as np @@ -201,3 +202,25 @@ def domain_plotter(title, mapextent, domains, background=None, resolution='1.0', alpha=0.7, transform=rotated_pole) return plt.gcf() + +def plot_cities(site_list): + ''' + Given a list of names of cities, plot them on the current axes + Must be included in the Natural Earth populated places + ''' + + # load populated places data from natural Earth + shpfilename = shpreader.natural_earth(resolution='10m', category='cultural', + name='populated_places') + reader = shpreader.Reader(shpfilename) + places = reader.records() + + # filter only sites in the list above + sites_shapes = [place for place in places if place.attributes['NAME'] in site_list ] + + ax = plt.gca() + + ax.scatter([point.attributes['LONGITUDE'] for point in sites_shapes], + [point.attributes['LATITUDE'] for point in sites_shapes], + c='red', + transform=ccrs.Geodetic()) From a34811729603073f211e3d3606b4055374e448b2 Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Mon, 16 Aug 2021 16:08:07 +0100 Subject: [PATCH 8/9] modify to be the new proposals for CNCP runs --- lib/catnip/test_plotter.py | 12 ++++++------ lib/catnip/um_domain.py | 4 ++-- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/lib/catnip/test_plotter.py b/lib/catnip/test_plotter.py index 5807dc6..5115f73 100644 --- a/lib/catnip/test_plotter.py +++ b/lib/catnip/test_plotter.py @@ -11,20 +11,20 @@ # 12lm domain-24 domain_12k = um_domain.UmDomain(delta_lat=0.11, delta_lon=0.11, frstlata=-24.0, - frstlona=328.0, + frstlona=310.0, polelata=pole_lat, polelona=pole_lon, global_row_length=700, global_rows=550) domain_4k_v2 = um_domain.UmDomain(delta_lat=0.0405, delta_lon=0.0405, - frstlata=0., frstlona=343.0, + frstlata=-2., frstlona=340.0, polelata=pole_lat, polelona=pole_lon, global_row_length=500, global_rows=432) # plot all of these on one map -mapextent = [80, 195, -10, 45] # x1, x2, y1, y2 +mapextent = [50, 178, -10, 45] # x1, x2, y1, y2 -title = '4 & 12 km domains. ' +title = ' New 4 & 12 km domains. ' fig = um_domain.domain_plotter(title, mapextent, [domain_4k_v2, domain_12k], background='topo', resolution='0.25', central_longitude=180) @@ -40,8 +40,8 @@ # and again for 4 km and inner sites site_list = ['Shanghai', 'Guangzhou', 'Shenzhen', 'Wuhan'] -mapextent = [108, 136, 18, 40] # x1, x2, y1, y2 -title = '4 km domain' +mapextent = [104, 132, 12, 36] # x1, x2, y1, y2 +title = 'New 4 km domain' fig = um_domain.domain_plotter(title, mapextent, [domain_4k_v2], background='topo', resolution='0.25', central_longitude=180) diff --git a/lib/catnip/um_domain.py b/lib/catnip/um_domain.py index 66115ca..4a35ff7 100644 --- a/lib/catnip/um_domain.py +++ b/lib/catnip/um_domain.py @@ -194,7 +194,7 @@ def domain_plotter(title, mapextent, domains, background=None, resolution='1.0', domain.minlat, domain.minlat] - p = ax.plot(x, y, marker='o', transform=rotated_pole) + p = ax.plot(x, y, marker=None, transform=rotated_pole) color = p[0].get_color() if plot_origin: # Add a marker for the origin of the domain @@ -206,7 +206,7 @@ def domain_plotter(title, mapextent, domains, background=None, resolution='1.0', def plot_cities(site_list): ''' Given a list of names of cities, plot them on the current axes - Must be included in the Natural Earth populated places + Must be included in the Natural Earth populated places ''' # load populated places data from natural Earth From 30c60f7d94c87ad9177feb50d090cb69650d7faa Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Mon, 16 Aug 2021 16:38:45 +0100 Subject: [PATCH 9/9] tweak to align it to original grids --- lib/catnip/{test_plotter.py => cncp_plots.py} | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) rename lib/catnip/{test_plotter.py => cncp_plots.py} (94%) diff --git a/lib/catnip/test_plotter.py b/lib/catnip/cncp_plots.py similarity index 94% rename from lib/catnip/test_plotter.py rename to lib/catnip/cncp_plots.py index 5115f73..a6e437f 100644 --- a/lib/catnip/test_plotter.py +++ b/lib/catnip/cncp_plots.py @@ -11,12 +11,12 @@ # 12lm domain-24 domain_12k = um_domain.UmDomain(delta_lat=0.11, delta_lon=0.11, frstlata=-24.0, - frstlona=310.0, + frstlona=309.96, polelata=pole_lat, polelona=pole_lon, global_row_length=700, global_rows=550) domain_4k_v2 = um_domain.UmDomain(delta_lat=0.0405, delta_lon=0.0405, - frstlata=-2., frstlona=340.0, + frstlata=-2.025, frstlona=340.003, polelata=pole_lat, polelona=pole_lon, global_row_length=500, global_rows=432)