Hi AndyH and other cfplotters:
I have precipitation on a rotated regional-climate-model grid, and I can get it to plot rather well with cfplot and cfview. But I can’t figure out how to avoid the bilinear interpolation. The original rotated grid has 34 by 27 values (ndim2 by ndmin1), but when I plot it, I would like to use nearest-neighbor sampling on the unrotated pixel coordinates. See the attached plots of what cfplot does, and what it looks like unrotated with nearest-neighbor sampling (with a different colour bar).
Can you help?
Patrick
cfplot
unrotated without cfplot and with a different colour bar
Hi Patrick,
I’m not quite sure what you mean about ‘avoid the bilinear interpolation’. Do you mean plot the data in rotated space and put a map on top? If so go into cfview → setup → map and select rotated as the map projection.
Cheers
Andy
Hi Andy:
I don’t really want to plot the data in rotated space. I want to plot the data (which is on a rotated grid) in lat/lon space. This I can do, like in the first plot I showed. But in that first plot, it is not plotting with nearest-neighbor sampling, but instead with bilinear interpolation. Ideally it would show a rotated and stretched regular checkerboard of 34 by 27 values. An unrotated, non-map-projected version of similar data is shown in the second plot in the last post, showing the 34 by 27 values with nearest-neighbor sampling and without bilinear interpolation.
I just tried to select ‘rotated’ in the cfview options for setup → map. Thanks for the suggestion. But that option isn’t there. I tried then to selected proj=‘rotated’ in cfplot.map(). That plotted the map without the lat/lon map projection. I don’t think that this is what I would like to do.
Maybe there’s some way to do the plot I suggest with the cfplot.bfill() routine? Or maybe there’s another way?
Patrick
Hi Patrick,
Could you let me know where the data is for these plots so I can take a look. I would imagine I need to do a modification in the bfill routine as you suggest.
Andy
Hi Andy:
Thanks. I put an example CORDEX subset file (with a rotated grid) on JASMIN in:
~pmcguire/for_andyh/
Patrick
Hi Patrick,
Does the following look to be what you’d like? I had to override the units and set them to ‘km’ as cf-python didn’t like the units of percentdiff when extracting the data from the field.
Yes, that is what I want.
How did you do it? Did you change your source code?
I used python code similar to this:
import cf
import cfplot as cfp
cfp.levs(min=-30, max=30, step=5)
cfp.mapset(lonmin=-80, lonmax=-70, latmin=-17, latmax=-7)
cfp.cscale('precip_diff_12lev')
xmargin = 0.03
rcp = 'rcp85'
cfp.gopen(user_position=**True**)
models=['MPI-M-MPI-ESM-LR','MOHC-HadGEM2-ES','NCC-NorESM1-M']
for i in range(0,3):
cfp.gpos(xmin=i*0.25 + xmargin,xmax=(i+1)*0.25 - xmargin, ymin=0.05, ymax=0.45)
f = cf.read('../CORDEX/pr_SAM-22_'+models[i]+'_'+rcp+'_r1i1p1_GERICS-REMO2015_v1_diff_2090-2099_2010-2019_timmean_percent.valleys+ancash2.nc')[0]
cfp.con(f, lines=False, titles=False, colorbar=**False**, title=models[i])
cfp.gclose()
Yes, I changed the cf-plot source code. I’ll need to do some more testing before I release it as a new version in case there’s any side effects for other plots. I’ll probably release it on Monday as I’m working from my laptop today and don’t have easy access to the release method.
Hi Patrick,
Just released the fixed code in cf-plot 3.2.15 to github and pypi. conda-forge will get the update automatically in a day or so. You have to use blockfill_fast=True option to cfp.con.
Hi Andy:
Excellent news! Thank you!
I will try it out as soon as I can.
Patrick
Hi AndyH:
I tried the rotated grid file with your new cfplot code with the same rotated-grid data file, and turning the interpolation off works fine (after overriding units). Thanks!
But with a 2nd file (without a time axis), the blockfill_fast=True option doesn’t work. It says "ValueError: Field.coordinate() can't return 2 constructs". See below.
Do I need to add a time axis somehow to get it to work? If so, do you know an easy way to do that? Or is there another way to get it to work without adding a time axis?
Patrick
import cf
import cfplot
datapath='~pmcguire/for_andyh/pr_SAM-22_MPI-M-MPI-ESM-LR_rcp26_r1i1p1_GERICS-REMO2015_v1_diff_2090-2099_2010-2019_timmean_percent.valleys+ancash2.nc'
f=cf.read(datapath)[0]
f.override_units('km', inplace=True)
f
<CF Field: precipitation_flux(time(1), grid_latitude(37), grid_longitude(34)) km>
#this works:
cfplot.con(f,blockfill_fast=True,lines=False)
datapath='~pmcguire/for_andyh/orog_SAM-22_MOHC-HadGEM2-ES_rcp85_r1i1p1_ICTP-RegCM4-7_v0_fx.valleys+[ancash2.nc](http://ancash2.nc/)'
f=cf.read(datapath)[0]
f
<CF Field: surface_altitude(projection_y_coordinate(39), projection_x_coordinate(37)) m>
#this doesn't work:
cfplot.con(f,blockfill_fast=True,lines=False)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/home/users/pmcguire/.local/lib/python3.10/site-packages/cfplot/cfplot.py", line 937, in con
if f.coord('X').has_bounds() and f.coord('Y').has_bounds():
File "/apps/jasmin/jaspy/miniconda_envs/jaspy3.10/m3-4.9.2/envs/jaspy3.10-m3-4.9.2-r20220721/lib/python3.10/site-packages/cf/mixin/fielddomain.py", line 3006, in coord
return self.coordinate(
File "/apps/jasmin/jaspy/miniconda_envs/jaspy3.10/m3-4.9.2/envs/jaspy3.10-m3-4.9.2-r20220721/lib/python3.10/site-packages/cf/mixin/fielddomain.py", line 1791, in coordinate
return self._construct(
File "/apps/jasmin/jaspy/miniconda_envs/jaspy3.10/m3-4.9.2/envs/jaspy3.10-m3-4.9.2-r20220721/lib/python3.10/site-packages/cf/mixin/fielddomain.py", line 200, in _construct
return self._default(
File "/apps/jasmin/jaspy/miniconda_envs/jaspy3.10/m3-4.9.2/envs/jaspy3.10-m3-4.9.2-r20220721/lib/python3.10/site-packages/cfdm/core/abstract/container.py", line 147, in _default
raise default
ValueError: Field.coordinate() can't return 2 constructs
Hi Patrick,
The time axis is probably not relevant here. The error is saying that there are two longitude and two latitude axes so cf-plot doesn’t know what to do as I’ve never seen data like this before.
list(f.coords())
[‘dimensioncoordinate0’,
‘dimensioncoordinate1’,
‘auxiliarycoordinate0’,
‘auxiliarycoordinate1’]
f.coord(‘dimensioncoordinate1’).X, f.coord(‘auxiliarycoordinate1’).X
(True, True)
f.coord(‘dimensioncoordinate0’).Y, f.coord(‘auxiliarycoordinate0’).Y
(True, True)
In your first dataset f.dump() gives:
Coordinate reference: grid_mapping_name:rotated_latitude_longitude
Coordinate conversion:grid_mapping_name = rotated_latitude_longitude
Coordinate conversion:grid_north_pole_latitude = 70.60000000000001
Coordinate conversion:grid_north_pole_longitude = -56.059999999999995
Coordinate conversion:north_pole_grid_longitude = 0.0
Dimension Coordinate: grid_longitude
Dimension Coordinate: grid_latitude
Auxiliary Coordinate: longitude
Auxiliary Coordinate: latitude
So this is what I have coded for in cf-plot.
In your second dataset f.dump() gives:
Coordinate reference: grid_mapping_name:oblique_mercator
Coordinate conversion:azimuth_of_central_line = 89.999999
Coordinate conversion:false_easting = -12500.0
Coordinate conversion:false_northing = -12500.0
Coordinate conversion:grid_mapping_name = oblique_mercator
Coordinate conversion:latitude_of_projection_origin = -22.0
Coordinate conversion:longitude_of_projection_origin = -59.0
Coordinate conversion:proj4_params = +proj=omerc +lat_0=-22.00 +alpha=89.999999 +lonc=-59.00 +x_0=-12500. +y_0=-12500. +ellps=sphere +a=6371229. +b=6371229. +units=m +no_defs
Coordinate conversion:scale_factor_at_projection_origin = 1.0
Datum:inverse_flattening = 0.0
Datum:semi_major_axis = 6371229.0
So this is something I haven’t seen before but contains various proj4 map conversion parameters. I have heard that Sadie is working on projections in cf-python so I will need to chat with her about this. There’s not much time left before the annual NCAS meeting so expect a delay in the response.
Hi AndyH
Thanks for looking at this.
Please do let me know what you find out from Sadie.
No rush.
Patrick
Hi Andy:
Have you been able to think a little more about this?
The rotated grid comes from CORDEX, which has local unrotated coordinates. That might be why there are two longitude and two latitude axes. Maybe it’s possible to tell cfplot to ignore the unrotated coordinates?
Patrick
Hi Andy:
I put a copy of the CMIP5 orography data in:
~/for_andyh/orog_fx_HadGEM2-ES_historical_r0i0p0.nc
It was copied from:
/badc/cmip5/data/cmip5/output1/MOHC/HadGEM2-ES/historical/fx/atmos/fx/r0i0p0/v20120215/orog/orog_fx_HadGEM2-ES_historical_r0i0p0.nc
This is using the settings below for South America west coast with your blockfill_fast=True code, but the ocean level doesn’t appear to be purple underflow (=0 elevation), like it is without blockfill_fast=True. Can you advise?
Thanks
Patrick
minlev = 0.0
steplev = 300.0
numlevs = 17 #predefined
maxlev = minlev+numlevs*steplev
cfp.mapset(lonmin=-74, lonmax=-65, latmin=-24, latmax=-16, proj='cyl')
cfp.levs(min=minlev, max=maxlev, step=steplev)
cfp.cscale('amwg_blueyellowred')
Hi Patrick, Yes you can do this as below.
f = cf.read('pmg2.nc')[0]
data = f.array
lons = f.coord('auxiliarycoordinate1').array
lats = f.coord('auxiliarycoordinate0').array
cfp.con(f=data, x=lons, y=lats, ptype=1)
Thanks, Andy!
That works without blockfill=true, but when I try it with blockfill=True, I get this error message:
>>> cfp.con(f=data, x=lons, y=lats, blockfill=True, ptype=1)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/home/users/pmcguire/.local/lib/python3.10/site-packages/cfplot/cfplot.py", line 952, in con
bfill(f=field_orig * fmult, x=x_orig, y=y_orig, clevs=clevs,
File "/home/users/pmcguire/.local/lib/python3.10/site-packages/cfplot/cfplot.py", line 4204, in bfill
xpts = np.append(xpts, x[ix] + (x[ix + 1] - x[ix]) / 2.0)
IndexError: index 39 is out of bounds for axis 0 with size 39
Any suggestions?
Patrick
Hi again, Andy:
Also, I am not sure which of the two issues you were addressing with the code sample from 1 hour ago. The purple underflow issue without blockfill=True? Or the “blockfill=True for the rotated grids” issue itself?
Patrick
Hi Andy:
I think your solution from about 2 hours ago was intended to get make the purple underflow appear for z<=0. Your code segment does appear to work for the South American region I indicated for the CORDEX rotated grid NETCDF without blockfill=True. Thanks! But I can’t get the purple underflow for z<=0 with the South American region I indicated with either:
1.the CMIP5 NETCDF with blockfill=True, or
2.the CORDEX rotated grid NETCDF with blockfill=True. This 2nd case of not working has the additional complication of index 39 being out of bounds for the axis of size 39.
Patrick
Hi Patrick,
You have two issues interleaving on the same ticket so It’s bound to get confusing!
Cordex data: The code I sent yesterday was to extract the longitudes and latitudes from the Cordex data and make a contour plot. I’m hoping to have a chat with Sadie later today on the unusual rotated coordinates and make some progress on that. The blockfill of the Cordex data doesn’t work at the moment in cf-plot.
Negative data plotting issue: I was thinking about this on the cycle home last night.
Let’s subspace the data to see what the values actually are:
mydata = f.subspace(X=cf.wi(-76, -72), Y=cf.wi(-24,-22))
mydata
<CF Field: surface_altitude(latitude(2), longitude(2)) m>
mydata.array
masked_array(
data=[[0., 0.],
[0., 0.]],
mask=False,
fill_value=1e+20,
dtype=float32)
So the values we have here are 0.0. In the colour scale below zero is coloured purple, 0.0 to 300.0 is dark blue and 300.0 to 600.0 is slightly less dark blue. Are we expecting negative numbers up to and including 0.0 to be coloured purple? Or are we expecting any values of 0.0 and above to 300.0 to be coloured dark blue? I would have thought that the latter was a more reasonable expectation. What do you think?
Of all the methods I think maybe the colour produced by blockfill_fast is the correct one.
