[ncl-talk] Contour plot with rotated grid shifted to the west
Karin Meier-Fleischer
meier-fleischer at dkrz.de
Tue Dec 1 08:08:31 MST 2015
Hi Guido,
I'm answering to ncl-talk so Mary and Dave can see my response to your
offline mail.
First, the script you send, had the wrong contour level min and max values.
Use the functions min and max to get more information about your data
values.
Also the colors you defined are not well chosen to start plotting your
data the first
time (most of the areas are filled with white color).
You have to use the values from your netCDF file for the rotated pole,
the given values in
the script are not correct.
Don't use /Res1 at tfDoNDCOverlay = True/ , but use /pv at lat2d = lat2d/ and
/pv at lon2d = lon2d/
instead. The png output looks ok, now.
I've cleaned the script so it is now easier to read. The script and the
plot are attached.
Bye,
Karin
-----------------
/Dear Karin, /
/I’m desperately trying to plot old data that I have in a rotated grid
on NCL. I stumbled upon your script that unrotated the coordinates and
plot the resulting data. However, it still doesn’t work for me, and I
don’t understand why, since Panoply is able to plot the same data./
//
/If you have some time please reply to this so that I can tell you more./
/Cheers//
///
/*
Guido Cioni*/
/*
*/
Am 01.12.15 um 10:33 schrieb Guido Cioni:
> I tried the script of Karin but still have issues. Since the bounds
> are not correctly computed the data are not correctly represented, as
> you can see in the following.
>
>
> I’m still convinced, though, that the data are correct since when I
> first ran the model I chose a rotation with center (-36 deg lon,45 deg
> lat) and left corner (-20 deg lon, -10 deg lat). Thus, the rotated
> pole has the correct coordinates (140 deg lon, 45 deg lat) since
> 140+36=180 and 45+45=90. If you take a look at rlat and rlon, that I
> printed with ncl_filedump, you can clearly see that in the rotated
> framework lat and lon both start from -10 (-9.962) and -20 degrees
> (360-340).
>
> Variable No. 0
> rlat:
>
> -9.962 -9.887 -9.812 -9.737 -9.662 -9.587
> -9.512 -9.437 -9.362 -9.287 -9.212 -9.137
> -9.062 -8.987 -8.912 -8.837 -8.762 -8.687
> -8.612 -8.537 -8.462 -8.387 -8.312 -8.237
> -8.162 -8.087 -8.012 -7.937 -7.862 -7.787
> -7.712 -7.637 -7.562 -7.487 -7.412 -7.337
> -7.262 -7.187 -7.112 -7.037 -6.962 -6.887
> -6.812 -6.737 -6.662 -6.587 -6.512 -6.437
> -6.362 -6.287 -6.212 -6.137 -6.062 -5.987
> -5.912 -5.837 -5.762 -5.687 -5.612 -5.537
> -5.462 -5.387 -5.312 -5.237 -5.162 -5.087
> -5.012 -4.937 -4.862 -4.787 -4.712 -4.637
> -4.562 -4.487 -4.412 -4.337 -4.262 -4.187
> -4.112 -4.037 -3.962 -3.887 -3.812 -3.737
> -3.662 -3.587 -3.512 -3.437 -3.362 -3.287
> -3.212 -3.137 -3.062 -2.987 -2.912 -2.837
> -2.762 -2.687 -2.612 -2.537 -2.462 -2.387
> -2.312 -2.237 -2.162 -2.087 -2.012 -1.937
> -1.862 -1.787 -1.712 -1.637 -1.562 -1.487
> -1.412 -1.337 -1.262 -1.187 -1.112 -1.037
> -0.9625 -0.8875 -0.8125 -0.7375 -0.6625 -0.5875
> -0.5125 -0.4375 -0.3625 -0.2875 -0.2125 -0.1375
> -0.6250E-01 0.1250E-01 0.8750E-01 0.1625 0.2375 0.3125
> 0.3875 0.4625 0.5375 0.6125 0.6875 0.7625
> 0.8375 0.9125 0.9875 1.063 1.138 1.213
> 1.288 1.363 1.438 1.513 1.588 1.663
> 1.738 1.813 1.888 1.963 2.038 2.113
> 2.188 2.263 2.338 2.413 2.488 2.563
> 2.638 2.713 2.788 2.863 2.938 3.013
> 3.088 3.163 3.238 3.313 3.388 3.463
> 3.538 3.613 3.688 3.763 3.838 3.913
> 3.988 4.063 4.138 4.213 4.288 4.363
> 4.438 4.513 4.588 4.663 4.738 4.813
> 4.888 4.963 5.038 5.113 5.188 5.263
> 5.338 5.413 5.488 5.563 5.638 5.713
> 5.788 5.863 5.938 6.013 6.088 6.163
> 6.238 6.313 6.388 6.463 6.538 6.613
> 6.688 6.763 6.838 6.913 6.988 7.063
> 7.138 7.213 7.288 7.363 7.438 7.513
> 7.588 7.663 7.738 7.813 7.888 7.963
> 8.038 8.113 8.188 8.263 8.338 8.413
> 8.488 8.563 8.638 8.713 8.788 8.863
> 8.938 9.013 9.088 9.163 9.238 9.313
> 9.388 9.463 9.538 9.613 9.688 9.763
> 9.838 9.913 9.988 10.06 10.14 10.21
> 10.29 10.36 10.44 10.51 10.59 10.66
> 10.74 10.81 10.89 10.96 11.04 11.11
> 11.19 11.26 11.34 11.41 11.49 11.56
> 11.64 11.71 11.79 11.86 11.94 12.01
> 12.09 12.16 12.24 12.31 12.39 12.46
> 12.54 12.61 12.69 12.76 12.84 12.91
> 12.99 13.06 13.14 13.21 13.29 13.36
> 13.44 13.51 13.59 13.66 13.74 13.81
> 13.89 13.96 14.04 14.11 14.19 14.26
> 14.34 14.41 14.49 14.56 14.64 14.71
> 14.79 14.86 14.94 15.01 15.09 15.16
> 15.24 15.31 15.39 15.46 15.54 15.61
> 15.69 15.76 15.84 15.91 15.99 16.06
> 16.14 16.21 16.29 16.36 16.44 16.51
> 16.59 16.66 16.74 16.81 16.89 16.96
> 17.04 17.11 17.19 17.26 17.34 17.41
> 17.49 17.56 17.64 17.71 17.79 17.86
> 17.94 18.01 18.09 18.16 18.24 18.31
> 18.39 18.46 18.54 18.61 18.69 18.76
> 18.84 18.91 18.99 19.06 19.14 19.21
> 19.29 19.36 19.44 19.51 19.59 19.66
> 19.74 19.81 19.89 19.96 20.04 20.11
> 20.19 20.26 20.34 20.41 20.49 20.56
> 20.64 20.71 20.79 20.86 20.94 21.01
> 21.09 21.16 21.24 21.31 21.39 21.46
> 21.54 21.61 21.69 21.76 21.84 21.91
> 21.99 22.06 22.14 22.21 22.29 22.36
> 22.44 22.51
>
> Variable No. 1
> rlon:
>
> 340.0 340.1 340.1 340.2 340.3 340.4
> 340.4 340.5 340.6 340.7 340.8 340.8
> 340.9 341.0 341.1 341.1 341.2 341.3
> 341.4 341.4 341.5 341.6 341.6 341.7
> 341.8 341.9 341.9 342.0 342.1 342.2
> 342.2 342.3 342.4 342.5 342.6 342.6
> 342.7 342.8 342.9 342.9 343.0 343.1
> 343.1 343.2 343.3 343.4 343.4 343.5
> 343.6 343.7 343.8 343.8 343.9 344.0
> 344.1 344.1 344.2 344.3 344.4 344.4
> 344.5 344.6 344.6 344.7 344.8 344.9
> 344.9 345.0 345.1 345.2 345.2 345.3
> 345.4 345.5 345.6 345.6 345.7 345.8
> 345.9 345.9 346.0 346.1 346.1 346.2
> 346.3 346.4 346.4 346.5 346.6 346.7
> 346.8 346.8 346.9 347.0 347.1 347.1
> 347.2 347.3 347.4 347.4 347.5 347.6
> 347.6 347.7 347.8 347.9 347.9 348.0
> 348.1 348.2 348.2 348.3 348.4 348.5
> 348.6 348.6 348.7 348.8 348.9 348.9
> 349.0 349.1 349.1 349.2 349.3 349.4
> 349.4 349.5 349.6 349.7 349.8 349.8
> 349.9 350.0 350.1 350.1 350.2 350.3
> 350.4 350.4 350.5 350.6 350.6 350.7
> 350.8 350.9 350.9 351.0 351.1 351.2
> 351.2 351.3 351.4 351.5 351.6 351.6
> 351.7 351.8 351.9 351.9 352.0 352.1
> 352.1 352.2 352.3 352.4 352.4 352.5
> 352.6 352.7 352.8 352.8 352.9 353.0
> 353.1 353.1 353.2 353.3 353.4 353.4
> 353.5 353.6 353.6 353.7 353.8 353.9
> 353.9 354.0 354.1 354.2 354.2 354.3
> 354.4 354.5 354.6 354.6 354.7 354.8
> 354.9 354.9 355.0 355.1 355.1 355.2
> 355.3 355.4 355.4 355.5 355.6 355.7
> 355.8 355.8 355.9 356.0 356.1 356.1
> 356.2 356.3 356.4 356.4 356.5 356.6
> 356.6 356.7 356.8 356.9 356.9 357.0
> 357.1 357.2 357.2 357.3 357.4 357.5
> 357.6 357.6 357.7 357.8 357.9 357.9
> 358.0 358.1 358.1 358.2 358.3 358.4
> 358.4 358.5 358.6 358.7 358.8 358.8
> 358.9 359.0 359.1 359.1 359.2 359.3
> 359.4 359.4 359.5 359.6 359.6 359.7
> 359.8 359.9 359.9 360.0 360.1 360.2
> 360.2 360.3 360.4 360.5 360.6 360.6
> 360.7 360.8 360.9 360.9 361.0 361.1
> 361.1 361.2 361.3 361.4 361.4 361.5
> 361.6 361.7 361.8 361.8 361.9 362.0
> 362.1 362.1 362.2 362.3 362.4 362.4
> 362.5 362.6 362.6 362.7 362.8 362.9
> 362.9 363.0 363.1 363.2 363.2 363.3
> 363.4 363.5 363.6 363.6 363.7 363.8
> 363.9 363.9 364.0 364.1 364.1 364.2
> 364.3 364.4 364.4 364.5 364.6 364.7
> 364.8 364.8 364.9 365.0 365.1 365.1
> 365.2 365.3 365.4 365.4 365.5 365.6
> 365.6 365.7 365.8 365.9 365.9 366.0
> 366.1 366.2 366.2 366.3 366.4 366.5
> 366.6 366.6 366.7 366.8 366.9 366.9
> 367.0 367.1 367.1 367.2 367.3 367.4
> 367.4 367.5 367.6 367.7 367.8 367.8
> 367.9 368.0 368.1 368.1 368.2 368.3
> 368.4 368.4 368.5 368.6 368.6 368.7
> 368.8 368.9 368.9 369.0 369.1 369.2
> 369.2 369.3 369.4 369.5 369.6 369.6
> 369.7 369.8 369.9 369.9 370.0 370.1
> 370.1 370.2 370.3 370.4 370.4 370.5
> 370.6 370.7 370.8 370.8 370.9 371.0
> 371.1 371.1 371.2 371.3 371.4 371.4
> 371.5 371.6 371.6 371.7 371.8 371.9
> 371.9 372.0 372.1 372.2 372.2 372.3
> 372.4 372.5 372.6 372.6 372.7 372.8
> 372.9 372.9 373.0 373.1 373.1 373.2
> 373.3 373.4 373.4 373.5 373.6 373.7
> 373.8 373.8 373.9 374.0 374.1 374.1
> 374.2 374.3 374.4 374.4 374.5 374.6
> 374.6 374.7 374.8 374.9 374.9 375.0
> 375.1 375.2 375.2 375.3 375.4 375.5
> 375.6 375.6 375.7 375.8 375.9 375.9
> 376.0 376.1 376.1 376.2 376.3 376.4
> 376.4 376.5 376.6 376.7 376.8 376.8
> 376.9 377.0 377.1 377.1 377.2 377.3
> 377.4 377.4 377.5 377.6 377.6 377.7
> 377.8 377.9 377.9 378.0 378.1 378.2
> 378.2 378.3 378.4 378.5 378.6 378.6
> 378.7 378.8 378.9 378.9 379.0 379.1
> 379.1 379.2 379.3 379.4 379.4 379.5
> 379.6 379.7 379.8 379.8 379.9 380.0
> 380.1 380.1 380.2 380.3 380.4 380.4
> 380.5 380.6 380.6 380.7 380.8 380.9
> 380.9 381.0 381.1 381.2 381.2 381.3
> 381.4 381.5 381.6 381.6 381.7 381.8
> 381.9 381.9 382.0 382.1 382.1 382.2
> 382.3 382.4 382.4 382.5 382.6 382.7
> 382.8 382.8 382.9 383.0 383.1 383.1
> 383.2 383.3 383.4 383.4 383.5 383.6
> 383.6 383.7 383.8 383.9 383.9 384.0
> 384.1 384.2 384.2 384.3 384.4 384.5
> 384.6 384.6 384.7 384.8 384.9 384.9
> 385.0 385.1 385.1 385.2 385.3 385.4
> 385.4 385.5 385.6 385.7 385.8 385.8
> 385.9 386.0 386.1 386.1 386.2 386.3
> 386.4 386.4 386.5 386.6 386.6 386.7
> 386.8 386.9 386.9 387.0 387.1 387.2
> 387.2 387.3 387.4 387.5 387.6 387.6
> 387.7 387.8 387.9 387.9 388.0 388.1
> 388.1 388.2 388.3 388.4 388.4 388.5
> 388.6 388.7 388.8 388.8 388.9 389.0
> 389.1 389.1 389.2 389.3 389.4 389.4
> 389.5 389.6 389.6 389.7 389.8 389.9
> 389.9 390.0 390.1 390.2 390.2 390.3
> 390.4 390.5 390.6 390.6 390.7 390.8
> 390.9 390.9 391.0 391.1 391.1 391.2
> 391.3 391.4 391.4 391.5 391.6 391.7
> 391.8 391.8 391.9 392.0 392.1 392.1
> 392.2 392.3 392.4 392.4 392.5 392.6
> 392.6 392.7 392.8 392.9
>
> However, the first image that I attached was referring to an older
> simulation, this is the correct one.
>
>
> Strangely enough, Panoply is able to plot the data just using
> rlon,rlat and the center of the rotated pole provided in the data file
> (I didn’t change anything to produce the following picture).
>
>
> Guido Cioni
> http://guidocioni.altervista.org
>
>
>
>
>> Il giorno 01 dic 2015, alle ore 09:52, Guido Cioni
>> <guidocioni at gmail.com <mailto:guidocioni at gmail.com>> ha scritto:
>>
>> Hey David,
>> yes, the image that I attached in the first email was supposed to
>> represents the final result. I know the “real” coordinate bounds just
>> because I wrote a binary file directly from the postprocessing of the
>> model with the real 2-D lat-lon coordinates in the native grid. Then,
>> I printed the minimum and maximum of those data. Funny thing is that
>> I wasn’t able to use these coordinates for my plot since they were
>> somewhat randomly organized, still don’t know why (probably some
>> fortran way of writing data…). However, the numbers that I attached
>> in the first plot were not correct. These are the correct ones:
>>
>> lon_min-78.64264
>> lon_max25.00719
>> lat_min27.55382
>> lat_max67.51249
>>
>> It puzzles me that they are not the same that you chose by
>> replicating the plot “by eye”. And also, you changed the center of
>> the map: I thought that the rotated pole coordinates that I had in
>> the file were the only correct data :)
>> Using the script that Karin wrote, I should be able to create
>> lat2d-lon2d array and unrotate the data, so that there is no need to
>> play around with the map projection, right?
>>
>>
>> Guido Cioni
>> http://guidocioni.altervista.org <http://guidocioni.altervista.org/>
>>
>>
>>
>>
>>> Il giorno 30 nov 2015, alle ore 23:54, David Brown <dbrown at ucar.edu
>>> <mailto:dbrown at ucar.edu>> ha scritto:
>>>
>>> Hi Guido,
>>> Here's my take on this so far. You enclosed an image labelled "Bolam
>>> Model". Is that plot supposed to represent the actual grid of your
>>> data? I am asking because it does not seem possible to reconcile the
>>> coordinate information that you have provided with this plot. I am
>>> attaching a couple of plots and a modified version of your script (it
>>> just draws the map, since we do not have the data). I can
>>> approximately duplicate the map projection of the Bolam plot by eye,
>>> but I have to use different numbers than those you have supplied. The
>>> script first draws a map based on the rotated lat/lon extent that you
>>> supplied, after conversion into "real" lat/lon values. Then it draws
>>> an approximation of the Bolam map.
>>>
>>> Mary supplied you with the NCL script from Karin that converts from
>>> rotated coordinates to regular lat/lon. I used a C code that does the
>>> same thing. You could use either method to create your own 2d lat/lon
>>> coordinates for this grid. If you would like the C code let me know.
>>>
>>> I am curious how you obtained the "real" corner values you used:
>>>
>>> Res1 at tfDoNDCOverlay=True
>>> Res1 at mpLimitMode = "Corners"
>>> Res1 at mpLeftCornerLatF = 30.35458
>>> Res1 at mpLeftCornerLonF = -75.14440
>>> Res1 at mpRightCornerLatF = 65.1125
>>> Res1 at mpRightCornerLonF = 13.87071
>>>
>>> These do not seem to me to be correct either.
>>>
>>> It would help if you listed the complete set of rotated lat/lon
>>> values. You could do it using ncl_filedump:
>>>
>>> ncl_filedump -v rlat,rlon complete.nc
>>>
>>> -dave
>>>
>>>
>>>
>>>
>>>
>>>
>>> On Mon, Nov 30, 2015 at 2:18 PM, Mary Haley <haley at ucar.edu
>>> <mailto:haley at ucar.edu>> wrote:
>>>> Correct, NCL can't automatically "unrotate" the coordinates for you
>>>> based on
>>>> the rotated coordinates.
>>>>
>>>> Karin Meier-Fleischer of DKRZ wrote a couple of functions that
>>>> "unrotate"
>>>> the coordinates for you, and I've taken this example and her data
>>>> file to
>>>> show how this might work for your data.
>>>>
>>>> See the attached script. The data file can be downloaded from:
>>>>
>>>> wget ftp.ucar.edu
>>>> <http://ftp.ucar.edu/>:/pub/scd/haley/tas_mon_1961-1990_rectilinear_grid_2D.nc
>>>>
>>>> The rotated grid actually doesn't look that "rotated", because the
>>>> center
>>>> lat/lon actually ends up being 0.0. But, hopefully you can use the
>>>> "unrot_lat" and "unrot_lon" functions to generate your own 2D lat/lon
>>>> arrays.
>>>>
>>>> --Mary
>>>>
>>>>
>>>> On Mon, Nov 30, 2015 at 1:31 PM, Guido Cioni <guidocioni at gmail.com
>>>> <mailto:guidocioni at gmail.com>> wrote:
>>>>>
>>>>> I don't have the original 2D lat/lon values, that's the problem.
>>>>> However,
>>>>> your script won't take into consideration the fact that
>>>>> coordinates are
>>>>> rotated, right? My rlon spans from approximately 340 to 392
>>>>> degrees, while
>>>>> rlat from - 10 to +22.. Of course these are not the real
>>>>> coordinates to be
>>>>> used for plotting the map, since they are expressed in the rotated
>>>>> framework. Do you see what I mean?
>>>>> I would like to provide you with the data file but is fairly big
>>>>> (gigs)
>>>>>
>>>>> Guido,
>>>>>
>>>>> I looked at your script further, and I think the issue is with the map
>>>>> corners you've selected. You are using the 1D rotated coordinate
>>>>> arrays to
>>>>> choose the corners, but this is likely not correct. What you
>>>>> really need are
>>>>> the original 2D lat/lon values, as the two examples illustrated,
>>>>> but I don't
>>>>> see that you have this information on the file.
>>>>>
>>>>> Did you try plotting the data using the coordinate arrays? I've
>>>>> created a
>>>>> script for you to try, although it's untested because I don't have
>>>>> your data
>>>>> file.
>>>>>
>>>>> Thanks,
>>>>>
>>>>> --Mary
>>>>>
>>>>>
>>>>> On Mon, Nov 30, 2015 at 11:45 AM, Guido Cioni
>>>>> <guidocioni at gmail.com <mailto:guidocioni at gmail.com>>
>>>>> wrote:
>>>>>>
>>>>>> Sure, here it is. I think that the problem is that ncl is
>>>>>> expecting 2-D
>>>>>> latitudes and longitudes while I have 1-D rlon and rlat..
>>>>>>
>>>>>> Yes I tried commenting the centers.
>>>>>>
>>>>>> m300382 at mlogin103% ncl_filedump complete.nc
>>>>>> Copyright (C) 1995-2015 - All Rights Reserved
>>>>>> University Corporation for Atmospheric Research
>>>>>> NCAR Command Language Version 6.3.0
>>>>>> The use of this software is governed by a License Agreement.
>>>>>> See http://www.ncl.ucar.edu/ for more details.
>>>>>>
>>>>>> Variable: f
>>>>>> Type: file
>>>>>> filename: complete
>>>>>> path: complete.nc
>>>>>> file global attributes:
>>>>>> CDI : Climate Data Interface version 1.4.7
>>>>>> (http://code.zmaw.de/projects/cdi)
>>>>>> Conventions : CF-1.0
>>>>>> history : Fri Nov 27 13:33:35 2015: cdo -f nc copy
>>>>>> 2014011700_000.grib2 2014011700_001.grib2 2014011700_002.grib2
>>>>>> 2014011700_003.grib2 2014011700_004.grib2 2014011700_005.grib2
>>>>>> 2014011700_006.grib2 2014011700_007.grib2 2014011700_008.grib2
>>>>>> 2014011700_009.grib2 2014011700_010.grib2 2014011700_011.grib2
>>>>>> 2014011700_012.grib2 2014011700_013.grib2 2014011700_014.grib2
>>>>>> 2014011700_015.grib2 2014011700_016.grib2 2014011700_017.grib2
>>>>>> 2014011700_018.grib2 2014011700_019.grib2 2014011700_020.grib2
>>>>>> 2014011700_021.grib2 2014011700_022.grib2 2014011700_023.grib2
>>>>>> 2014011700_024.grib2 2014011700_025.grib2 2014011700_026.grib2
>>>>>> 2014011700_027.grib2 2014011700_028.grib2 2014011700_029.grib2
>>>>>> 2014011700_030.grib2 2014011700_031.grib2 2014011700_032.grib2
>>>>>> 2014011700_033.grib2 2014011700_034.grib2 2014011700_035.grib2
>>>>>> 2014011700_036.grib2 2014011700_037.grib2 2014011700_038.grib2
>>>>>> 2014011700_039.grib2 2014011700_040.grib2 2014011700_041.grib2
>>>>>> 2014011700_042.grib2 2014011700_043.grib2 2014011700_044.grib2
>>>>>> 2014011700_045.grib2 2014011700_046.grib2 2014011700_047.grib2
>>>>>> CDO : Climate Data Operators version 1.4.7
>>>>>> (http://code.zmaw.de/projects/cdo)
>>>>>> dimensions:
>>>>>> ncl_scalar = 1
>>>>>> rlon = 706
>>>>>> rlat = 434
>>>>>> lev = 20
>>>>>> time = 73 // unlimited
>>>>>> variables:
>>>>>> double rlon ( rlon )
>>>>>> long_name : longitude in rotated pole grid
>>>>>> units : degrees
>>>>>> standard_name : grid_longitude
>>>>>> axis : X
>>>>>>
>>>>>> double rlat ( rlat )
>>>>>> long_name : latitude in rotated pole grid
>>>>>> units : degrees
>>>>>> standard_name : grid_latitude
>>>>>> axis : Y
>>>>>>
>>>>>> character rotated_pole ( ncl_scalar )
>>>>>> grid_mapping_name : rotated_latitude_longitude
>>>>>> grid_north_pole_latitude : 44.999996
>>>>>> grid_north_pole_longitude : 144
>>>>>>
>>>>>> double lev ( lev )
>>>>>> long_name : pressure
>>>>>> units : Pa
>>>>>> axis : Z
>>>>>>
>>>>>> double time ( time )
>>>>>> units : hours since 2014-01-17 00:00:00
>>>>>> calendar : proleptic_gregorian
>>>>>>
>>>>>> float gh ( time, lev, rlat, rlon )
>>>>>> long_name : Geopotential Height
>>>>>> units : gpm
>>>>>> code : 5
>>>>>> grid_mapping : rotated_pole
>>>>>>
>>>>>> float t ( time, lev, rlat, rlon )
>>>>>> long_name : Temperature
>>>>>> units : K
>>>>>> grid_mapping : rotated_pole
>>>>>>
>>>>>> float u ( time, lev, rlat, rlon )
>>>>>> long_name : U velocity
>>>>>> units : m s**-1
>>>>>> code : 2
>>>>>> grid_mapping : rotated_pole
>>>>>>
>>>>>> float v ( time, lev, rlat, rlon )
>>>>>> long_name : V velocity
>>>>>> units : m s**-1
>>>>>> code : 3
>>>>>> grid_mapping : rotated_pole
>>>>>>
>>>>>> float w ( time, lev, rlat, rlon )
>>>>>> long_name : Vertical velocity
>>>>>> units : Pa s**-1
>>>>>> code : 8
>>>>>> grid_mapping : rotated_pole
>>>>>>
>>>>>> float pv ( time, lev, rlat, rlon )
>>>>>> long_name : Potential vorticity
>>>>>> units : K m**2 kg**-1 s**-1
>>>>>> code : 14
>>>>>> grid_mapping : rotated_pole
>>>>>>
>>>>>> float mterh ( time, rlat, rlon )
>>>>>> long_name : Model terrain height
>>>>>> units : m
>>>>>> code : 7
>>>>>> grid_mapping : rotated_pole
>>>>>>
>>>>>> Il 30 nov 2015 6:28 PM, "Mary Haley" <haley at ucar.edu
>>>>>> <mailto:haley at ucar.edu>> ha scritto:
>>>>>>>
>>>>>>> Guido,
>>>>>>>
>>>>>>> I'm hoping Dave will weigh in here, but it would help if we could at
>>>>>>> least see the output from:
>>>>>>>
>>>>>>> ncl_filedump pv.nc
>>>>>>>
>>>>>>> Sometimes the global attributes of the file will give some valuable
>>>>>>> information about the correct map projection parameters.
>>>>>>>
>>>>>>> Did you actually try commenting out the mpCenterLatF and
>>>>>>> mpCenterLonF
>>>>>>> settings? In the script you included, these two resources are
>>>>>>> still being
>>>>>>> set.
>>>>>>>
>>>>>>> It's possible that the data has already been rotated for you.
>>>>>>>
>>>>>>> --Mary
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> On Mon, Nov 30, 2015 at 1:33 AM, Guido Cioni
>>>>>>> <guidocioni at gmail.com <mailto:guidocioni at gmail.com>>
>>>>>>> wrote:
>>>>>>>>
>>>>>>>> Hi all,
>>>>>>>> I’m trying to produce a contour plot with data on a rotated
>>>>>>>> grid. The
>>>>>>>> original data look like that
>>>>>>>>
>>>>>>>>
>>>>>>>> Since I have the rotated north pole latitude (44.999996) and
>>>>>>>> longitude
>>>>>>>> (144.0) in my data file I used the second example of the
>>>>>>>> support page on
>>>>>>>> rotated grid.
>>>>>>>>
>>>>>>>> Res1 at tfDoNDCOverlay=True
>>>>>>>> Res1 at mpLimitMode = "Corners"
>>>>>>>> Res1 at mpLeftCornerLatF = lat(0)
>>>>>>>> Res1 at mpLeftCornerLonF = lon(0)
>>>>>>>> Res1 at mpRightCornerLatF = lat(nlat-1)
>>>>>>>> Res1 at mpRightCornerLonF = lon(nlon-1)
>>>>>>>> Res1 at mpCenterLatF = 90 - north_pole_lat
>>>>>>>> Res1 at mpCenterLonF = 180 + north_pole_lon
>>>>>>>>
>>>>>>>> However I get this result ( I believe this is Australia):
>>>>>>>>
>>>>>>>>
>>>>>>>> After several attempts I was able to plot the data using as
>>>>>>>> bounds the
>>>>>>>> “real” coordinates of the rotated projection, without
>>>>>>>> specifying the center,
>>>>>>>> e.g. :
>>>>>>>>
>>>>>>>> Res1 at tfDoNDCOverlay=True
>>>>>>>> Res1 at mpLimitMode = "Corners"
>>>>>>>> Res1 at mpLeftCornerLatF = 30.35458
>>>>>>>> Res1 at mpLeftCornerLonF = -75.14440
>>>>>>>> Res1 at mpRightCornerLatF = 65.1125
>>>>>>>> Res1 at mpRightCornerLonF = 13.87071
>>>>>>>>
>>>>>>>> But the contour plot is still slightly shifted to the west as
>>>>>>>> you can
>>>>>>>> see in this image:
>>>>>>>>
>>>>>>>>
>>>>>>>> I know that in the data file my coordinates span from
>>>>>>>> rlon={340,392.875} and rlat={-9.962499,22.512501}.
>>>>>>>>
>>>>>>>> What am I missing?
>>>>>>>> Is there any additional documentation on rotated grid?
>>>>>>>>
>>>>>>>> Attached is my NCL script (data are too big).
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> Guido Cioni
>>>>>>>> http://guidocioni.altervista.org
>>>>>>>> <http://guidocioni.altervista.org/>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> ncl-talk mailing list
>>>>>>>> ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu>
>>>>>>>> List instructions, subscriber options, unsubscribe:
>>>>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>>>>>
>>>>>>>
>>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> ncl-talk mailing list
>>>> ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu>
>>>> List instructions, subscriber options, unsubscribe:
>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>
>>> <pv_test.000001.png><pv_test.000002.png><plot_pv-mod.ncl>
>>
>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
--
Dipl. Geophys. Karin Meier-Fleischer
Visualization, NCL
Application Support
Deutsches Klimarechenzentrum GmbH (DKRZ)
Bundesstrasse 45a - D20146 Hamburg - Germany
Phone: +49 (0)40 460094 126
Fax: +49 (0)40 460094 270
E-Mail: meier-fleischer at dkrz.de
URL: www.dkrz.de
Geschäftsführer: Prof. Dr. Thomas Ludwig
Sitz der Gesellschaft: Hamburg
Amtsgericht Hamburg HRB 39784
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151201/03f5c7ee/attachment-0001.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/png
Size: 409496 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151201/03f5c7ee/attachment-0004.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/png
Size: 152982 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151201/03f5c7ee/attachment-0005.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/png
Size: 244755 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151201/03f5c7ee/attachment-0006.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pv_test.png
Type: image/png
Size: 109443 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151201/03f5c7ee/attachment-0007.png
-------------- next part --------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;------------------------------------------------------------
;-- set global constants
;------------------------------------------------------------
pi = 4.0*atan(1.)
deg2rad = pi/180.
rad2deg = 45./atan(1.)
fillval = -99999.9
;------------------------------------------------------------
;-- Function: create_lon2d(y,x)
;-- Description: create 2d lon array from 1D rotlon, rotlat
;------------------------------------------------------------
undef("create_lon2d")
function create_lon2d(rotlat[*]:numeric, rotlon[*]:numeric)
local x, i
begin
x = new((/dimsizes(rotlat),dimsizes(rotlon)/),typeof(rotlat))
; print("Function create_lon2d: dimsizes of x: "+dimsizes(x))
do i=0,dimsizes(rotlat)-1
x(i,:) = rotlon
end do
return(x)
end
;------------------------------------------------------------
;-- Function: create_lat2d(y,x)
;-- Description: create 2d lat array from 1D rotlon,rotlat
;------------------------------------------------------------
undef("create_lat2d")
function create_lat2d(rotlat[*]:numeric, rotlon[*]:numeric)
local y, i
begin
y = new((/dimsizes(rotlat),dimsizes(rotlon)/),typeof(rotlat))
; print("Function create_lat2d: dimsizes of y: "+dimsizes(y))
do i=0,dimsizes(rotlon)-1
y(:,i) = rotlat
end do
return(y)
end
;------------------------------------------------------------
;-- Function: unrot_lon(rotlat,rotlon,pollat,pollon)
;-- Description: transform rotated longitude to longitude
;------------------------------------------------------------
undef("unrot_lon")
function unrot_lon( rotlat:numeric, rotlon:numeric, pollat[1]:numeric, pollon[1]:numeric )
local rotlat, rotlon, pollon, pollat, lon, s1, c1, s2, c2, rlo, rla, i, tmp1, tmp2
begin
lon = fillval
lon at _FillValue = fillval
if (any(dimsizes(rotlat) .ne. dimsizes(rotlon)) .and. \
(dimsizes(dimsizes(rotlat)).ne.1 .or. dimsizes(dimsizes(rotlon)).ne.1)) then
print("Function unrot_lon: unrot_lon: rotlat and rotlon dimensions do not match")
return(lon)
end if
if (dimsizes(dimsizes(rotlat)).eq.1 .and. dimsizes(dimsizes(rotlon)).eq.1) then
rla = create_lat2d(rotlat,rotlon) ;-- create 2D latitude array
rlo = create_lon2d(rotlat,rotlon) ;-- create 2D longitude array
else
rla = rotlat
rlo = rotlon
end if
rla = rla*deg2rad ;-- convert from degree to radians
rlo = rlo*deg2rad ;-- convert from degree to radians
lon := (/rlo/) ;-- reassign lon
lon at _FillValue=fillval
s1 = sin(pollat*deg2rad)
c1 = cos(pollat*deg2rad)
s2 = sin(pollon*deg2rad)
c2 = cos(pollon*deg2rad)
tmp1 = s2*(-s1*cos(rlo)*cos(rla)+c1*sin(rla))-c2*sin(rlo)*cos(rla)
tmp2 = c2*(-s1*cos(rlo)*cos(rla)+c1*sin(rla))+s2*sin(rlo)*cos(rla)
lon = atan(tmp1/tmp2)*rad2deg
; lon = where(lon.lt.0, lon+360,lon)
; lon = where(lon.gt.180, lon-180,lon)
lon at units = "degrees_east"
print("Function unrot_lon: min/max "+sprintf("%8.4f", min(lon(0,:)))+" "+sprintf("%8.4f", max(lon(0,:))))
delete([/rlo,rlo,c1,s1,c2,s2,tmp1,tmp2/])
return(lon)
end
;------------------------------------------------------------
;-- Function: unrot_lat(rotlat,rotlon,pollat,pollon)
;-- Description: transform rotated latitude to latitude
;------------------------------------------------------------
undef("unrot_lat")
function unrot_lat( rotlat:numeric, rotlon:numeric, pollat[1]:numeric, pollon[1]:numeric )
local rotlat, rotlon, pollon, pollat, lat, s1, c1, rlo, rla, i
begin
lat = fillval
lat at _FillValue = fillval
if (any(dimsizes(rotlat) .ne. dimsizes(rotlon)) .and. \
(dimsizes(dimsizes(rotlat)).ne.1 .or. dimsizes(dimsizes(rotlon)).ne.1)) then
print("Function unrot_lat: rotlat and rotlon dimensions do not match")
return(lat)
end if
if (dimsizes(dimsizes(rotlat)).eq.1 .and. dimsizes(dimsizes(rotlon)).eq.1) then
rla = create_lat2d(rotlat,rotlon) ;-- create 2D latitude array
rlo = create_lon2d(rotlat,rotlon) ;-- create 2D longitude array
else
rla = rotlat
rlo = rotlon
end if
rla = rla*deg2rad ;-- convert from degree to radians
rlo = rlo*deg2rad ;-- convert from degree to radians
lat := (/rla/) ;-- reassign lat
lat at _FillValue=fillval
s1 = sin(pollat*deg2rad)
c1 = cos(pollat*deg2rad)
lat = s1*sin(rla)+c1*cos(rla)*cos(rlo)
lat = asin(lat)*rad2deg
lat at units = "degrees_north"
print("Function unrot_lat: min/max "+sprintf("%8.4f", min(lat(:,0)))+" "+sprintf("%8.4f", max(lat(:,0))))
delete([/rlo,rla,c1,s1/])
return(lat)
end
;----------------
;-- MAIN
;----------------
begin
; Open model level output file
; File = addfile ("/Volumes/mistral_scratch/bolam/complete.nc","r")
File = addfile("/scratch/m/m300382/bolam/complete.nc","r")
rlat = File->rlat
rlon = File->rlon
rotpole = File->rotated_pole
pollat = rotpole at grid_north_pole_latitude
pollon = rotpole at grid_north_pole_longitude
; read data
pv=File->pv
; delete(pv at units)
; delete(pv at long_name)
lon2d = unrot_lon(rlat, rlon, pollat, pollon)
lat2d = unrot_lat(rlat, rlon, pollat, pollon)
;-- set 2D lat/lons
pv at lat2d := lat2d ;-- 2D latitudes
pv at lon2d := lon2d ;-- 2D longitudes
;-- calculate the min and max lat/lons for the map plot
minlat := min(lat2d) ;-- retrieve minimum latitude value
minlon := min(lon2d) ;-- retrieve maximum latitude value
maxlat := max(lat2d) ;-- retrieve minimum longitude value
maxlon := max(lon2d) ;-- retrieve maximum longitude value
print("----> MinLatF: "+minlat+" MaxLatF: "+maxlat)
print("----> MinLonF: "+minlon+" MaxLonF: "+maxlon)
nlat = dimsizes(lat2d(:,0))
nlon = dimsizes(lon2d(0,:))
; create plot
wks = gsn_open_wks("png","pv_test")
Res1 = True
Res1 at gsnDraw = False
Res1 at gsnFrame = False
Res1 at gsnAddCyclic = False ;-- data are not global, don't add lon cyclic point
Res1 at pmTickMarkDisplayMode = "Always" ;-- draw tickmarks
Res1 at mpDataBaseVersion = "MediumRes" ;-- use finer database
Res1 at mpMinLatF = minlat - 1. ;-- set min lat
Res1 at mpMaxLatF = maxlat + 1. ;-- set max lat
Res1 at mpMinLonF = minlon - 1. ;-- set min lon
Res1 at mpMaxLonF = maxlon + 1. ;-- set max lon
Res1 at lbLabelBarOn = False ; turn off individual cb's
Res1 at cnFillOn = True ; do color fill
Res1 at cnLinesOn = False
Res1 at cnInfoLabelOn = False
Res1 at cnLineLabelsOn = False
Res1 at cnFillPalette = "rainbow"
; colors=(/"white","lightyellow","darkorange","red"/)
; rgb = span_named_colors(colors,False)
; gsn_define_colormap(wks,rgb)
plot=new(2, graphic)
Res1 at tiMainString="17/01/2014 21 UTC"
plot(0) = gsn_csm_contour_map(wks,pv({22},{35000},:,:),Res1)
Res1 at tiMainString="18/01/2014 10 UTC"
plot(1) = gsn_csm_contour_map(wks,pv({43},{35000},:,:),Res1)
;-- create panel
pnlres = True
pnlres at txString = ""
pnlres at gsnPanelLabelBar = True
pnlres at gsnMaximize = True
gsn_panel(wks,plot,(/2,1/),pnlres)
end
More information about the ncl-talk
mailing list