[ncl-talk] Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates
Dave Allured - NOAA Affiliate
dave.allured at noaa.gov
Fri Mar 5 15:56:31 MST 2021
I checked the data source. This problem goes all the way back to the data
set that is published on the CEDA Archive website. Their 12 km subset over
UK contains the same invalid coordinates in grid_latitude and
grid_longitude variables. I asked them to fix their data set, so let's see
what happens there.
Note that there is no problem with Iris; they are just copying the invalid
coordinates into the file that SBR made.
Thanks to Karin for the regridding solution.
On Fri, Mar 5, 2021 at 3:29 AM Karin Meier-Fleischer via ncl-talk <
ncl-talk at mailman.ucar.edu> wrote:
> Hi,
>
> I found the solution to regrid OSGB pre-projected data to lonlat grid
> with CDO.
>
> To remap the data to a lonlat grid over GB, first create a grid file
> using the topo operator:
>
> cdo -f nc -sellonlatbox,-14,6,48,62 -remapbil,r1440x720 -topo topo.nc
>
> Get the grid information using the griddes operator. You have to modify
> the grid information file (see uploaded grid_OSGB.txt file): delete the
> upper curvilinear part, add 'nvertex = 2', and the line
> proj_params = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717
> +x_0=400000 +y_0=-100000 +ellps=airy
> +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m
> +no_defs"
>
> Then remap your data file to the lonlat grid:
>
> cdo -remapbil,topo.nc -setgrid,grid_OSGB.txt -selname,tas \
> tas_rcp85_land-rcm_uk_12km.nc outfile.nc
>
> It is important to use the setgrid oprator, too.
>
> -Karin
>
> Am 01.03.21 um 13:08 schrieb S Br:
> > Hi Karin,
> > Thanks for trying this. The data is certainly working fine with Python.
> > Please find the script here.
> > I don't have expertise in Python but took help from one of my colleagues.
> >
> > import iris
> > import iris.analysis
> > import iris.plot as iplt
> > import matplotlib.pyplot as plt
> > import iris.quickplot as qplt
> >
> > file1 = 'tasmax_rcp85_land-rcm_uk_12km_time10.nc
> > <http://tasmax_rcp85_land-rcm_uk_12km_time10.nc>'
> > cubes1 = iris.load(file1)
> > print(cubes1)
> > #tas1 = cubes1[1]
> > #print(tas1)
> >
> > plt.figure(figsize=(10,10))
> > plt.subplot(1, 2, 1)
> > iplt.pcolormesh(cubes1[1][0,0,:,:])
> > plt.title('Original grid')
> > ax = plt.gca()
> > ax.coastlines()
> > ax.gridlines()
> > plt.tight_layout()
> > plt.show()
> > ax = plt.gca()
> > ax.coastlines(resolution='50m')
> > ax.gridlines()
> > plt.tight_layout()
> > plt.show()
> >
> > Thanks.
> >
> > On Mon, Mar 1, 2021 at 11:12 AM Karin Meier-Fleischer via ncl-talk
> > <ncl-talk at mailman.ucar.edu <mailto:ncl-talk at mailman.ucar.edu>> wrote:
> >
> > As Dave and Dennis already wrote there is something wrong with the
> > grid_latitude and grid_longitude values. I was not able to get the
> > correct plot neither with NCL, Python nor CDO.
> >
> > How did you create a correct plot with Python? Did you use Iris? It
> > would be very interesting to see your Python script. Is it possible
> to
> > get it?
> >
> > -Karin
> >
> >
> > Am 26.02.21 um 17:38 schrieb Dave Allured - NOAA Affiliate via
> ncl-talk:
> > > Progress. This proves that arrays grid_latitude and
> > grid_longitude are
> > > wrong in the input file. This is no longer an NCL problem,
> > except for
> > > possibly trying to reproject the data using projection
> > parameters. I
> > > think Dennis already said this projection is not supported in
> > NCL, so
> > > that would be difficult.
> > >
> > > I recommend asking Iris support for help to generate the correct
> > lat and
> > > lon arrays. You can show them the invalid coordinates directly,
> > leaving
> > > NCL out of the loop, by sampling the output from:
> > >
> > > ncdump -v grid_latitude,grid_longitude file.nc
> > <http://file.nc> <http://file.nc <http://file.nc>>
> > >
> > > It may also be possible to get the correct coordinate arrays by
> some
> > > other method. You might ask the original data provider if they
> can
> > > supply the correct lat and lon arrays to match their data. If
> > you can
> > > get corrected arrays, it would be easy to retrofit them into your
> > > temperature data file.
> > >
> > >
> > > On Fri, Feb 26, 2021 at 2:34 AM S Br <sbr.climate at gmail.com
> > <mailto:sbr.climate at gmail.com>
> > > <mailto:sbr.climate at gmail.com <mailto:sbr.climate at gmail.com>>>
> wrote:
> > >
> > > Hi,
> > > The min and max of f->grid_latitude and f->grid_longitude
> prints
> > > wrong information. You can see that the bounding box is not
> > > surrounding the UK.
> > >
> > > (0) min=-0.8710618968191547 max=12.94984143216418
> > > (0) min=-18.24145632530888 max=-6.580251255557555
> > >
> > > Thanks.
> > >
> > > On Thu, Feb 25, 2021 at 9:24 PM Dave Allured - NOAA Affiliate
> > > <dave.allured at noaa.gov <mailto:dave.allured at noaa.gov>
> > <mailto:dave.allured at noaa.gov <mailto:dave.allured at noaa.gov>>>
> wrote:
> > >
> > > Okay. We need to evaluate whether the lat and lon
> > coordinates
> > > in your file are correct. I see Dennis's other reply,
> > but I am
> > > suggesting a different approach using only attached
> > coordinates,
> > > and nothing else. If the Iris package generated your
> file
> > > correctly, then nothing else should be needed.
> > >
> > > Please show the min and max of f->grid_latitude and
> > > f->grid_longitude. Don't show the whole array. The
> command
> > > *printMinMax (var, 0)* is useful for this. If these
> > coordinate
> > > sets are correct, then the min and max should be a
> > bounding box
> > > surrounding the UK.
> > >
> > > NCL works best with supplied coordinate arrays, which is
> the
> > > method I am recommending here. Python is doing something
> > > different, perhaps recalculating based on the projection
> > > information. I don't really care about that, I just want
> to
> > > help you get the normal NCL method working.
> > >
> > >
> > > On Thu, Feb 25, 2021 at 1:20 PM S Br
> > <sbr.climate at gmail.com <mailto:sbr.climate at gmail.com>
> > > <mailto:sbr.climate at gmail.com
> > <mailto:sbr.climate at gmail.com>>> wrote:
> > >
> > > Hi Dave,
> > > Please see the attached plot, this is how it looks
> > like if
> > > we use the above statements.
> > > It obviously doesn't plot over the UK.
> > >
> > >
> > > On Thu, Feb 25, 2021 at 6:59 PM Dave Allured - NOAA
> > > Affiliate <dave.allured at noaa.gov
> > <mailto:dave.allured at noaa.gov>
> > > <mailto:dave.allured at noaa.gov
> > <mailto:dave.allured at noaa.gov>>> wrote:
> > >
> > > I see what's going on. From your other e-mail,
> > the data
> > > file contains two different sets of coordinates.
> > There
> > > are the projection coordinates in meters. Those
> > are not
> > > latitudes and longitudes, therefore not
> > sufficient for
> > > NCL. But the geographic coordinates are also
> > provided.
> > > Double coordinates are frequently encountered in
> GIS
> > > related applications. Try this:
> > >
> > > tas at lat2d=f->grid_latitude
> > > tas at lon2d=f->grid_longitude
> > >
> > >
> > > On Thu, Feb 25, 2021 at 11:42 AM S Br
> > > <sbr.climate at gmail.com
> > <mailto:sbr.climate at gmail.com> <mailto:sbr.climate at gmail.com
> > <mailto:sbr.climate at gmail.com>>>
> > > wrote:
> > >
> > > If I print tas at lat2d and tas at lon2d, they
> give the
> > > following information.
> > >
> > > Variable: taslat
> > > Type: double
> > > Total Size: 896 bytes
> > > 112 values
> > > Number of Dimensions: 1
> > > Dimensions and sizes: [112]
> > > Coordinates:
> > > (0) -102000
> > > (1) -90000
> > > (2) -78000
> > > (3) -66000
> > > (4) -54000
> > > (5) -42000
> > > (6) -30000
> > > (7) -18000
> > > (8) -6000
> > > (9) 6000
> > > (10) 18000
> > > (11) 30000
> > > (12) 42000
> > > (13) 54000
> > > (14) 66000
> > > (15) 78000
> > > (16) 90000
> > > (17) 102000
> > > (18) 114000
> > > (19) 126000
> > > (20) 138000
> > > (21) 150000
> > > (22) 162000
> > > (23) 174000
> > > (24) 186000
> > > (25) 198000
> > > (26) 210000
> > > (27) 222000
> > > (28) 234000
> > > (29) 246000
> > > (30) 258000
> > > (31) 270000
> > > (32) 282000
> > > (33) 294000
> > > (34) 306000
> > > (35) 318000
> > > (36) 330000
> > > (37) 342000
> > > (38) 354000
> > > (39) 366000
> > > (40) 378000
> > > (41) 390000
> > > (42) 402000
> > > (43) 414000
> > > (44) 426000
> > > (45) 438000
> > > (46) 450000
> > > (47) 462000
> > > (48) 474000
> > > (49) 486000
> > > (50) 498000
> > > (51) 510000
> > > (52) 522000
> > > (53) 534000
> > > (54) 546000
> > > (55) 558000
> > > (56) 570000
> > > (57) 582000
> > > (58) 594000
> > > (59) 606000
> > > (60) 618000
> > > (61) 630000
> > > (62) 642000
> > > (63) 654000
> > > (64) 666000
> > > (65) 678000
> > > (66) 690000
> > > (67) 702000
> > > (68) 714000
> > > (69) 726000
> > > (70) 738000
> > > (71) 750000
> > > (72) 762000
> > > (73) 774000
> > > (74) 786000
> > > (75) 798000
> > > (76) 810000
> > > (77) 822000
> > > (78) 834000
> > > (79) 846000
> > > (80) 858000
> > > (81) 870000
> > > (82) 882000
> > > (83) 894000
> > > (84) 906000
> > > (85) 918000
> > > (86) 930000
> > > (87) 942000
> > > (88) 954000
> > > (89) 966000
> > > (90) 978000
> > > (91) 990000
> > > (92) 1002000
> > > (93) 1014000
> > > (94) 1026000
> > > (95) 1038000
> > > (96) 1050000
> > > (97) 1062000
> > > (98) 1074000
> > > (99) 1086000
> > > (100) 1098000
> > > (101) 1110000
> > > (102) 1122000
> > > (103) 1134000
> > > (104) 1146000
> > > (105) 1158000
> > > (106) 1170000
> > > (107) 1182000
> > > (108) 1194000
> > > (109) 1206000
> > > (110) 1218000
> > > (111) 1230000
> > >
> > >
> > > Variable: taslon
> > > Type: double
> > > Total Size: 656 bytes
> > > 82 values
> > > Number of Dimensions: 1
> > > Dimensions and sizes: [82]
> > > Coordinates:
> > > (0) -210000
> > > (1) -198000
> > > (2) -186000
> > > (3) -174000
> > > (4) -162000
> > > (5) -150000
> > > (6) -138000
> > > (7) -126000
> > > (8) -114000
> > > (9) -102000
> > > (10) -90000
> > > (11) -78000
> > > (12) -66000
> > > (13) -54000
> > > (14) -42000
> > > (15) -30000
> > > (16) -18000
> > > (17) -6000
> > > (18) 6000
> > > (19) 18000
> > > (20) 30000
> > > (21) 42000
> > > (22) 54000
> > > (23) 66000
> > > (24) 78000
> > > (25) 90000
> > > (26) 102000
> > > (27) 114000
> > > (28) 126000
> > > (29) 138000
> > > (30) 150000
> > > (31) 162000
> > > (32) 174000
> > > (33) 186000
> > > (34) 198000
> > > (35) 210000
> > > (36) 222000
> > > (37) 234000
> > > (38) 246000
> > > (39) 258000
> > > (40) 270000
> > > (41) 282000
> > > (42) 294000
> > > (43) 306000
> > > (44) 318000
> > > (45) 330000
> > > (46) 342000
> > > (47) 354000
> > > (48) 366000
> > > (49) 378000
> > > (50) 390000
> > > (51) 402000
> > > (52) 414000
> > > (53) 426000
> > > (54) 438000
> > > (55) 450000
> > > (56) 462000
> > > (57) 474000
> > > (58) 486000
> > > (59) 498000
> > > (60) 510000
> > > (61) 522000
> > > (62) 534000
> > > (63) 546000
> > > (64) 558000
> > > (65) 570000
> > > (66) 582000
> > > (67) 594000
> > > (68) 606000
> > > (69) 618000
> > > (70) 630000
> > > (71) 642000
> > > (72) 654000
> > > (73) 666000
> > > (74) 678000
> > > (75) 690000
> > > (76) 702000
> > > (77) 714000
> > > (78) 726000
> > > (79) 738000
> > > (80) 750000
> > > (81) 762000
> > >
> > > On Thu, Feb 25, 2021 at 6:26 PM Dave Allured
> > - NOAA
> > > Affiliate <dave.allured at noaa.gov
> > <mailto:dave.allured at noaa.gov>
> > > <mailto:dave.allured at noaa.gov
> > <mailto:dave.allured at noaa.gov>>> wrote:
> > >
> > > It is starting to look like there is an
> > error in
> > > the attached coordinates. What are the
> > min and
> > > max values of tas at lat2d and tas at lon2d?
> > >
> > >
> > > On Thu, Feb 25, 2021 at 10:45 AM S Br
> > > <sbr.climate at gmail.com
> > <mailto:sbr.climate at gmail.com>
> > > <mailto:sbr.climate at gmail.com
> > <mailto:sbr.climate at gmail.com>>> wrote:
> > >
> > > Hi Dennis,
> > > The data draws on a wrong location
> > and the
> > > map boundary is not aligning.
> > > It should plot over the UK which is
> about
> > > 50-60N.
> > >
> > > Thanks.
> > >
> > > On Thu, Feb 25, 2021 at 5:41 PM
> > Dennis Shea
> > > <shea at ucar.edu <mailto:shea at ucar.edu>
> > <mailto:shea at ucar.edu <mailto:shea at ucar.edu>>> wrote:
> > >
> > > It helps to give NCL the correct
> > lat/lon
> > > array information.
> > >
> > > See attached
> > >
> > > On Thu, Feb 25, 2021 at 10:23 AM
> Dave
> > > Allured - NOAA Affiliate via
> ncl-talk
> > > <ncl-talk at mailman.ucar.edu
> > <mailto:ncl-talk at mailman.ucar.edu>
> > > <mailto:ncl-talk at mailman.ucar.edu
> > <mailto:ncl-talk at mailman.ucar.edu>>> wrote:
> > >
> > > That script looks pretty
> > good. You
> > > can see from the error
> > messages that
> > > NCL is trying to
> > understand @lat2d
> > > and @lon2d, but something is
> not
> > > quite right with them. We
> > need to
> > > see printVarSummary for tas,
> > > tas at lat2d, and tas at lon2d.
> It is
> > > possible that dimensions are
> > > mismatched between data and
> > > coordinates, for some reason.
> > >
> > > There is also this note in
> > that doc
> > > page: "If your data variable
> > is a
> > > two-dimensional (2D) array
> > ordered
> > > lon x lat, you will need to
> > reorder
> > > it to be lat x lon." Please
> > check
> > > and see if this is the
> problem.
> > >
> > >
> > > On Thu, Feb 25, 2021 at 10:00
> > AM S
> > > Br <sbr.climate at gmail.com
> > <mailto:sbr.climate at gmail.com>
> > > <mailto:sbr.climate at gmail.com
> > <mailto:sbr.climate at gmail.com>>> wrote:
> > >
> > >
> > > Hi Dave,
> > > Many thanks for your
> reply.
> > > Sorry I missed to provide
> the
> > > details.
> > > I have run the following
> NCL
> > > script and came across
> the
> > > errors given below.
> > > You can see the attached
> plot
> > > that it draws value
> > throughout
> > > the globe. It is supposed
> > to be
> > > drawn only over the UK.
> > > Am I doing
> something wrong?
> > >
> > > ;********
> > >
> > srcFileName="tas_rcp85_land-rcm_uk_12km_01_mon_198012-208011.nc
> > <http://tas_rcp85_land-rcm_uk_12km_01_mon_198012-208011.nc>
> > >
> > <http://tas_rcp85_land-rcm_uk_12km_01_mon_198012-208011.nc
> > <http://tas_rcp85_land-rcm_uk_12km_01_mon_198012-208011.nc>>"
> > >
> > > ;---Read File
> > > f=addfile(srcFileName,"r")
> > > print(f)
> > >
> > > ;---Get the source
> > lat/lon grid
> > > tas=f->tas(0,0,:,:)
> > >
> > tas at lat2d=f->projection_y_coordinate
> > >
> > tas at lon2d=f->projection_x_coordinate
> > >
> > > ;---Plot data
> > > wks =
> > gsn_open_wks("x11","map")
> > >
> > > res = True
> > > res at gsnMaximize =
> > True ;
> > > maximize plot in frame
> > > res at cnFillOn =
> > True ;
> > > turn on contour fill
> > > res at cnLinesOn =
> > False ;
> > > turn off contour fill
> > > res at gsnAddCyclic =
> > False
> > > ; res at mpMinLatF =
> 48.0
> > > ; only plot 30S
> > to 30N
> > > ; res at mpMaxLatF =
> 60.0
> > > ; res at mpMinLonF =
> -15.0
> > > ; only plot 30S to
> 30N
> > > ; res at mpMaxLonF =
> 6.0
> > > res at mpDataBaseVersion =
> > > "MediumRes"
> > > res at mpLandFillColor =
> > "white"
> > >
> > > plot =
> > >
> > gsn_csm_contour_map(wks,tas(:,:),res)
> > >
> > > ********************
> > > (0)
> is_valid_latlon2d_attr:
> > > Warning: The 'lat2d'
> > attribute
> > > must either be
> > > (0) the same dimension
> > sizes as
> > > the data, or one element
> > larger
> > > in both directions.
> > > (0) Your data will most
> > likely
> > > not be overlaid on the map
> > > correctly.
> > > (0) check_for_y_lat_coord:
> > > Warning: Data either does
> not
> > > contain
> > > (0) a valid latitude
> > coordinate
> > > array or doesn't contain
> > one at all.
> > > (0) A valid latitude
> > coordinate
> > > array should have a
> 'units'
> > > (0) attribute equal to
> one of
> > > the following values:
> > > (0) 'degrees_north'
> > > 'degrees-north'
> > 'degree_north'
> > > 'degrees north'
> 'degrees_N'
> > > 'Degrees_north' 'degree_N'
> > > 'degreeN' 'degreesN' 'deg
> > north'
> > > (0)
> is_valid_latlon2d_attr:
> > > Warning: The 'lon2d'
> > attribute
> > > must either be
> > > (0) the same dimension
> > sizes as
> > > the data, or one element
> > larger
> > > in both directions.
> > > (0) Your data will most
> > likely
> > > not be overlaid on the map
> > > correctly.
> > > (0) check_for_lon_coord:
> > > Warning: Data either does
> not
> > > contain
> > > (0) a valid longitude
> > coordinate
> > > array or doesn't contain
> > one at all.
> > > (0) A valid longitude
> > coordinate
> > > array should have a
> 'units'
> > > (0) attribute equal to
> one of
> > > the following values:
> > > (0) 'degrees_east'
> > > 'degrees-east'
> 'degree_east'
> > > 'degrees east' 'degrees_E'
> > > 'Degrees_east' 'degree_E'
> > > 'degreeE' 'degreesE' 'deg
> > east'
> > >
> > > On Thu, Feb 25, 2021 at
> > 4:17 PM
> > > Dave Allured - NOAA
> Affiliate
> > > <dave.allured at noaa.gov
> > <mailto:dave.allured at noaa.gov>
> > >
> > <mailto:dave.allured at noaa.gov <mailto:dave.allured at noaa.gov>>>
> > > wrote:
> > >
> > > Does your file
> > contain its
> > > own lat and lon
> > coordinate
> > > arrays? You will
> need to
> > > check this yourself
> with
> > > ncdump or NCL plus
> print
> > > commands.
> > >
> > > If the file has its
> own
> > > coordinate arrays,
> > then do
> > > not regrid or
> reproject.
> > > Doing that will lose
> > detail
> > > and reduce the
> quality of
> > > the plotted image.
> > NCL is
> > > quite capable of
> directly
> > > displaying any data
> with
> > > attached coordinate
> > arrays,
> > > and it does not need
> > to know
> > > anything about the
> > > projection of the
> data.
> > > Please see
> > documentation for
> > > "Plotting data on a
> > map" for
> > > much more information
> > about
> > > aligning data on a
> > map plot.
> > >
> > > https://www.ncl.ucar.edu/Applications/plot_data_on_map.shtml
> > <https://www.ncl.ucar.edu/Applications/plot_data_on_map.shtml>
> > >
> > <https://www.ncl.ucar.edu/Applications/plot_data_on_map.shtml
> > <https://www.ncl.ucar.edu/Applications/plot_data_on_map.shtml>>
> > >
> > > If the coordinate
> > arrays in
> > > the file are two
> > > dimensional, then see
> > > example 3,
> "Curvilinear
> > > grid". When you have
> 2-D
> > > coordinates, they
> must be
> > > passed to the NCL plot
> > > routine with special
> > > attributes *lat2d* and
> > > *lon2d*, or
> alternatively
> > > *sfXArray* and
> > > *sfYArray* attached
> > to the
> > > plot resource
> variable.
> > > Note that 1-D
> > coordinates do
> > > not need this special
> > > treatment because
> > they are
> > > normally attached
> > directly
> > > to the data variable
> > to plot.
> > >
> > >
> > > On Thu, Feb 25, 2021
> > at 8:24
> > > AM S Br via ncl-talk
> > >
> > <ncl-talk at mailman.ucar.edu <mailto:ncl-talk at mailman.ucar.edu>
> > >
> > <mailto:ncl-talk at mailman.ucar.edu <mailto:
> ncl-talk at mailman.ucar.edu>>>
> > > wrote:
> > >
> > > Hi All,
> > > I have a set of
> > regional
> > > gridded data in
> the
> > > British National
> Grid
> > > (OSGB) spatial
> > > coordinate
> > system. I am
> > > looking to use
> > this data
> > > in NCL but the
> > > projection is not
> > >
> correctly displayed?
> > > Could it be
> > possible to
> > > regrid this data
> to
> > >
> Latitude-longitude in
> > > rotated pole
> > coordinates
> > > or to a regular
> > lat/lon
> > > grid.
> > >
> > > I would like to
> > mention
> > > that these data
> > sets are
> > > generated using
> > the Iris
> > > package in
> > Python. The
> > > data is correctly
> > > displayed in
> > Python but
> > > not in NCL.
> > >
> > > For your
> > convenience I
> > > have attached
> here a
> > > sample NetCDF
> file.
> > >
> > > Thank you.
> > > SBR
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20210305/71da52ba/attachment.html>
More information about the ncl-talk
mailing list