[ncl-talk] Regridding from OSGB spatial coordinate to Latitude-longitude in rotated pole coordinates
Karin Meier-Fleischer
meier-fleischer at dkrz.de
Fri Feb 26 01:42:30 MST 2021
Hi,
here, I weigh in to make some comments about CDO and the Python Iris
projections. Your question about the OSGB data in the CDO forum had me
running into the same problems days ago.
CDO usually can handle pre-projected grids using a gridfile with the
grid_mapping informaion, but in your case there is something missing to
create a correct gridfile. Transverse Mercator projection is quite special.
The python package Iris provides an extra class for the projection type
OSGB.
class iris.coord_systems.OSGBA
Specific transverse mercator projection on a specific ellipsoid
I'll see how Iris has solved the problem and get back to you, but maybe
not today;).
Ciao,
Karin
Am 25.02.21 um 22:07 schrieb Dennis Shea via ncl-talk:
> I don't understand:
>
> "This data *can't also read* correctly with Climate Data Operator (CDO)"
>
> To me something is not correct with the file if the CDO and NCL are
> having issues.
> Python is rendering a correct plot? Somebody has programmed around
> the issue.
> ===
>
> You can go back and forth between the lat/lon gridded coordinates and
> x/y as indicated by the equations at:
>
> *https://idlecoding.com/transverse-mercator-with-python/comment-page-1/#:~:text=%20Transverse%20Mercator%20With%20Python%20%201%20The,any%20code%20to%20transform%20an%20actual...%20More%20*
> <https://idlecoding.com/transverse-mercator-with-python/comment-page-1/#:~:text=%20Transverse%20Mercator%20With%20Python%20%201%20The,any%20code%20to%20transform%20an%20actual...%20More%20>
>
> But to me, these are both on the file.
>
> double grid_longitude(projection_y_coordinate,
> projection_x_coordinate) ;
> double grid_latitude(projection_y_coordinate,
> projection_x_coordinate) ;
>
> double projection_x_coordinate(projection_x_coordinate) ;
> double projection_y_coordinate(projection_y_coordinate) ;
>
> The information at:
>
> int transverse_mercator ;
> transverse_mercator:grid_mapping_name =
> "transverse_mercator" ;
> transverse_mercator:longitude_of_prime_meridian = 0. ;
> transverse_mercator:semi_major_axis = 6377563.396 ;
> transverse_mercator:semi_minor_axis = 6356256.909 ;
> transverse_mercator:longitude_of_central_meridian = -2. ;
> transverse_mercator:latitude_of_projection_origin = 49. ;
> transverse_mercator:false_easting = 400000. ;
> transverse_mercator:false_northing = -100000. ;
> transverse_mercator:scale_factor_at_central_meridian =
> 0.9996012717 ;
>
> is possibly useful, I have tried to use two of the attributes but it did
> not work.
>
>
>
>
> On Thu, Feb 25, 2021 at 1:08 PM S Br <sbr.climate at gmail.com
> <mailto:sbr.climate at gmail.com>> wrote:
>
> Hi Dennis,
> The plot is perfect with Python. Please see the attached image.
> The data I am talking about is generated using Python.
>
> This data can't also read correctly with Climate Data Operator (CDO)
>
> Thanks.
>
>
> On Thu, Feb 25, 2021 at 7:26 PM Dennis Shea <shea at ucar.edu
> <mailto:shea at ucar.edu>> wrote:
>
> Obviously, not over the UK but ......
>
> I think more information from you might be necessary.
> "transverse_mercator" is not supported by NCL options.
>
> You state: "The data is correctly displayed in Python but not in
> NCL."
>
> Well, what does the plot look like via Python?
>
> On Thu, Feb 25, 2021 at 12:19 PM Dennis Shea <shea at ucar.edu
> <mailto:shea at ucar.edu>> wrote:
>
> I **speculate** that the attributes associated with the
> "transverse_mercator" variable must be used.
> ====================================================
> ;---Get the source lat/lon grid
> tas = f->tasmax(0,0,:,:)
>
> lat2d = f->grid_latitude
> lon2d = f->grid_longitude
> printMinMax(lat2d,0)
> printMinMax(lon2d,0)
>
> itm = f->transverse_mercator
>
> print(itm)
> print("itm at latitude_of_projection_origin="+itm at latitude_of_projection_origin)
> print("itm at longitude_of_central_meridian="+itm at longitude_of_central_meridian)
>
> ; clarity
>
> lat_shift = itm at latitude_of_projection_origin
> lon_shift = itm at longitude_of_central_meridian
>
> ; apply 'shift'
>
> lat2d = lat2d + lat_shift
> lon2d = lon2d + lon_shift
> printMinMax(lat2d,0)
> printMinMax(lon2d,0)
>
> tas at lat2d = lat2d ; use 'shifted' coordinates
> tas at lon2d = lon2d
> ================
>
> See attached
>
>
> On Thu, Feb 25, 2021 at 12:00 PM Dave Allured - NOAA
> Affiliate via ncl-talk <ncl-talk at mailman.ucar.edu
> <mailto:ncl-talk at mailman.ucar.edu>> 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>>
> 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>> 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>> 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>> 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>> 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>> 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>"
>
> ;---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>>
> 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>
>
> 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>>
> 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
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at mailman.ucar.edu <mailto:ncl-talk at mailman.ucar.edu>
> List instructions, subscriber options, unsubscribe:
> https://mailman.ucar.edu/mailman/listinfo/ncl-talk
> <https://mailman.ucar.edu/mailman/listinfo/ncl-talk>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at mailman.ucar.edu
> List instructions, subscriber options, unsubscribe:
> https://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
More information about the ncl-talk
mailing list