Friday, May 3, 2013

Tutorial: Overlay image on basemap image using GrADS

This tutorial teaches you how to save GrADS output onto an already existing image, a basemap.  This allows you to bring more color and vibrance to your weather map by adding, for example, topographic maps underneath of your weather data.  The use of basemaps can transform the image on the left to the image on the right.
Before Basemap
After Basemap

Before we get started with the tutorial, it is recommended that you learn how to make your own basemaps (see tutorial on making basemaps here).  However, if you want to skip learning how to make basemaps for now, I will provide the basemap needed for this tutorial.  The basemap provided will be the same basemap used in the above example.  To download this basemap to your computer, simply right click on the image below and choose "Save Image As..." and save it to your GrADS folder.

Save this image, as your basemap for this tutorial


So, now that you have your basemap, lets get started.  Plotting data onto a basemap is very simple, and requires very little extra code.  All you need to do is add a few extra options onto the 'printim' command, when you are outputting your data to an image.
  • Normal: 'printim img.png png x800 y600'
  • Basemap: 'printim img.png png x800 y600 -b basemap.png -t 0'
The options -b and -t, respectively, stand for 'basemap' and 'transparency'.  After the -b you must indicate the path to the basemap image, and after the -t you must specify which color to make transparent.  You can find out this number by querying the shading information through the 'q shades' command, or simply specifying this number through the 'set ccols' and 'set clevs' command.

And that's pretty much all there is to overlaying your data onto a basemap.

There are however a few things to watch out for.  Since all you are essentially doing when you include the -b option, is printing your GrADS plot onto a preexisting image, you need to make sure all the image dimensions match.  To do this there are really only three things that you need to remember:
  1. Your GrADS page area (parea) matches the page area on your basemap
  2. Your image dimensions are equal (you can get your image dimensions by right clicking your image and looking at 'properties')
  3. Your lat/lon boundaries must be the same between the basemap domain and the GrADS domain
If these three conditions are not satisfied, you will end up with an overlay that does not coincide with your basemap (see below image for an example of what this looks like)

Incorrect!

Now, that you have an idea on how to do it, here are the correct dimensions for the basemap provided here.

Image Dimensions: x=800px, y=600px
Parea: '0.5 10.0 0.5 7.5'
latitude: 22, 51
longitude: -128, -65

So if you were to use these dimensions in a script to plot radar, your code might look like:

     'set mpdset hires'
     'set lat 22 51'

     'set lon -128 -65'                                                        ;*Depending on the dataset, you may need to add 360 to this
     'set display color white'
     'clear'

     'set parea 0.5 10.0 0.5 7.5'                                                  ;*Note the page area is equal to the basemap
     'set gxout shaded'
     'set ccols 0 14 4 11 5 13 19 7 12 8 2 6'                                ;*Note that the 0 is the specified transparent color
     'set clevs  5 10 15 20 25 30 35 40 45 50 55 60 65'
     'd radar'
     'printim img.png png -b basemap.png -t 0 x800 y600'
       ;*x800, y800 correspond to the image dimensions


That about does it for this tutorial, you should now be able to overlay images onto preexisting basemaps.  And just for fun, and so you know, your background image does not necessarily have to be a map.


I hope you enjoyed, and found this tutorial useful!




15 comments:

  1. Hi,

    I'm plotting my snow maps above a background image of europe, but when there is no snow being calculated by the model, my background image turns white (so the color background isn't visible). Is there a special code for getting back the color map, even when there is no snow being calculated?

    Thanks in advance

    ReplyDelete
  2. Two things to check, make sure that you set your contour levels ('set clevs') and your colors ('set ccols') Also, instead of setting the lower level of snow to zero try just above zero (like 0.01). Then the last thing is to make sure that when you use the 'printim' command, you specify your transparent color (-t) to the color marked at your lowest contour level.

    Example:
    'set clevs 0.01 1 2 4 5'
    'set ccols 0 2 3 4 5'

    'printim img.png png -b basemap.png -t 0'

    This way, any snow value lower than 0.01 will be assigned a color with the numeric value 0, then when you plot the image on your basemap, any colors with the numeric value 0 will be masked out.

    If this does not help, you can always use the maskout function, which only plots values where the snow amount exceeds some prescribed level. See Here:

    http://gradsaddict.blogspot.com/2013/05/tutorialtwo-methods-to-plot-more-than.html

    Hope that helps!

    ReplyDelete
  3. Hello,

    I downloaded the script provided above and changed the last line in the script to reflect where I wanted to place the image output, the name of the file and the directory to where the basemap/overlay image was located respectively (after the -b).

    'printim /home/wxprob/Desktop/image.png png -b /home/wxprob/Desktop/basic.png -t 0 x800 y600'

    Unfortunately, I receive an error that says, 'PRINTIM error: background image open error.'

    What is it that I am doing wrong? This is on a linux operating system, if that helps.

    ReplyDelete
  4. The problem is simple, GrADS cannot find your background image, so double check your file path to the background image. If you are using a linux system, you may need more in front of the /home/ to get the full file path. Go to the folder with your image and type 'pwd' into the terminal, this should give you your full file path.

    ReplyDelete
  5. Hello from Italy, Fiumicino Airport
    Thanks a lot for all scripts over here. can you tell me if you know how to calculate reflectivity for Europe from Grads variables ? Thanks a lot Mario

    ReplyDelete
    Replies
    1. So, If you have the reflectivity data, either from a radar or model output, this isn't hard. There are a few tutorials here that can help you out. If you want to "calculate" a reflectivity estimate given other data, then I am afraid I can't help.

      Sorry I couldn't be of more help, I am unfamiliar with how to access data in Europe.

      Delete
    2. Hello and thanks for the answer. European data are GFS data. Since GFS don't have "reflectivity" unlike WRF, I see someone putting togherer "rain, ice and othere parameters"and create a virtual reflectivity. I just cannot find the sources. Mario

      Delete
  6. Greetings from Sub-Saharan!!

    I would like to make an ASCII file similar to lpoly_US.asc which I can use to mask out areas that are not within East Africa. How can I make such a file?

    ReplyDelete
  7. Hi,

    Is there anyone that has an opoly_hires.asc and lpoly_hires.asc for africa?

    Regards

    ReplyDelete
  8. Hello! I'm learning grads and wanted to let you know how awesome this site of yours is, it has taught me a ton of stuff already.

    I have one question... I am trying to plot grads data on a basemap using your colormaps script. I tried picking colormaps colors (for example 's3pcpn') that began with white and then put -t 0 in the printim command like shown above, but all I get is the shading with no map background whatsoever. Any advice?

    ReplyDelete
    Replies
    1. Hi Charlie, so I am pretty sure I know what your issue is but without the code in front of me I can't be 100%. So I'll give you a quick answer here, and then if that doesn't work send me the code and I can look at it more in depth.

      So the way I set up colormaps.gs is similar to how color.gs is set up in that the assigned colors don't overwrite the 16 default colors that GrADS uses. So, the first color in any map is numbered as "16", therefore the white color in the s3pcpn is number 16, not 0. My guess is that is you change the -t 0 to -t 16, you should fix the problem.

      Good luck!
      -Ted

      Delete
    2. Unfortunately, that didn't work, but I appreciate the suggestion. I was able to get it to work with -t 1 (probably because I have 'set display color white') for the standard ccols/clevs method but wasn't able to get it to work with the color.gs or colormaps.gs scripts. I have my code here:

      Basemap: https://github.com/WeatherTogether/model-charts/blob/master/usbasemap.gs

      Surface CAPE (the thing I'm having trouble with): https://github.com/WeatherTogether/model-charts/blob/master/sfccape.gs

      Thank you so much for your help!
      Charlie

      Delete
    3. Okay, I looked at your code and I have one other thing you can try. Try setting the minimum value of CAPE in the colormaps call to be 5 or some other low value instead of zero. What I think is happening is the code sets everything BELOW the minimum value to 16, and everything equal to the minimum value to 17. Since the minimum value of CAPE is zero. If that doesn't work I'll dig around a bit more.

      Delete
    4. Unfortunately, that didn't work either. One note I forgot to add: I can still see axes labels from my basemap with color.gs or colormaps.gs enabled (with -t 1) but can't see the topography on the basemap.

      Delete
    5. I figured it out! There were two white levels, one at 1 and one at 16. Thanks for your help, and thanks for writing such an easy-to-understand and useful blog.

      Delete