Monday, December 10, 2012

Raster calculations with GDAL and numpy: calculating GFS wid speed

Performing raster calculations is a frequent need. GDAL python bindings come with a utility program,  gdal_calc.py, that does that, but many times is useful to code them.
In this example, we will be able to calculate the wind speed field from the GFS model, from the given u and v speed.
 The calculated wind field. America is easily seen at the right side of the image

The data

In this example, we will use the GFS global model data, which is available at their web site. Just navigate to a date that pleases you and download a .grb file.
Files are in grib format, with a little problem when trying to view them in QGis: Coordinates are from 0 to 360º in longitude, while the usual thing is to put them from -180º to 180º. Everything will work without changing this, but you can do it after taking a look to this post.

Inside the grib file, there are many bands (hundreds). You can see their definition using the gdal command gdalifo, or the file with the extension *inv at the same place you downloaded the grib file. There are many wind fields, that correspond to different altitudes, all them indicated as u (east) and v components (north). Just take note of the band number of both at the desired level.

NumPy

The best of using the python gdal bindings, compared from those of the other programs, is that converting a raster file to a NumPy matrix is now very easy. 
NumPy is a python extension that allows to perform matrix calculations in a really easy and efficient way. It should be installed in your computer if GDAL python is installed.
Once you have two matrices A and B in NumPy, adding element by element to a C matrix would be as easy as:
C = A + B
and multiplying:
C = A * B
Note that is not a matrix multiplication (more info here), but this is what we usually need when working with raster images.

The code

The wind speed is the modulus of the two components, so 
spd = sqrt(u^2+v^2)
The "traditional way" to calculate it would be opening two bands and loop all the elements to create the new band. Instead, this can be done:
from osgeo import gdal
from osgeo.gdalnumeric import *
from osgeo.gdalconst import *

fileName = "gfs.tiff"
bandNum1 = 199
bandNum2 = 200

outFile = "out.tiff"

#Open the dataset
ds1 = gdal.Open(fileName, GA_ReadOnly )
band1 = ds1.GetRasterBand(bandNum1)
band2 = ds1.GetRasterBand(bandNum2)

#Read the data into numpy arrays
data1 = BandReadAsArray(band1)
data2 = BandReadAsArray(band2)

#The actual calculation
dataOut = numpy.sqrt(data1*data1+data2*data2)

#Write the out file
driver = gdal.GetDriverByName("GTiff")
dsOut = driver.Create("out.tiff", ds1.RasterXSize, ds1.RasterYSize, 1, band1.DataType)
CopyDatasetInfo(ds1,dsOut)
bandOut=dsOut.GetRasterBand(1)
BandWriteArray(bandOut, dataOut)

#Close the datasets
band1 = None
band2 = None
ds1 = None
bandOut = None
dsOut = None

Note that:
  1. gdalnumeric is used, because it's the easiest way to dump a raster into a NumPy matrix using BandReadAsArray and the inverse function with BandWriteArray.
  2. In this case, we are reading the bands 199 and 200, but you should check which bands are the ones to read in the file you download.
  3. CopyDatasetInfo will copy all the input file metadata into the output file. This includes the geotransform, the projection, the GCPs (there aren't in this case) and the other metadata.
  4. Don't forget to set the bands and datasets to None at the end to free them. In longer scripts this may be critical.

gdal_calc.py

If you want to use the command line function to do the same, you can try with:
gdal_calc.py -A gfs.tiff --A_band 199 -B gfs.tiff --B_band 200 --outfile out.tiff --calc="sqrt(A*A+B*B)"
It gives the same result.

43 comments:

  1. Hi,

    thanks for this post, very clear.

    I am a math teacher and, to introduce the concepts of vector field / ordinary differential equations,
    I'd like to plot a wind field over our french Brittany region. Here are the coordinates using the basemap module
    notations :

    urcrnrlon=-2.6, urcrnrlat=49.0,
    llcrnrlon=-5.34, llcrnrlat=47

    Yet I am a newbie in the geosciences area. So I think (tell me please if I am wrong) the workflow is the following :

    1. Get a grb file as you said,
    2. Extract from this file the geographic zone I am interested in (Brittany coasts)
    3. Generate the raster for the wind field (corresponding the gfs.tiff file you used).
    4. Make computations as you did in your script (essentially compute the modulus of the vector field)

    I have no idea with 2 and 3 ... (since 1 and 4 are trivial thanks to your post !)

    Last question : how does the gdalifo command works ?

    Thanks a lot !

    Pr. Yassine Patel
    yaspat.github.io/blog

    ReplyDelete
    Replies
    1. Hi Yassine,

      I think that using the Basemap library will make your life easier. This example will plot you the results in the frame defined when you create the Basemap instance: https://basemaptutorial.readthedocs.org/en/latest/plotting_data.html#contourf
      Of course, you have to change the data array to whatever you need (the module, in your case).
      I like using basemap, since you can show the wind direction too:
      https://basemaptutorial.readthedocs.org/en/latest/plotting_data.html#barbs
      or, even prettier:
      https://basemaptutorial.readthedocs.org/en/latest/plotting_data.html#quiver

      I'm happy to know that the blog posts can be used to teach!

      Delete
  2. Hi,
    can we create mosaic dataset in gdal or numpy?
    thanks

    ReplyDelete
  3. if gdal cannot read raster as array, how could we use raster calculator using gdal?

    ReplyDelete
    Replies
    1. Which data format are you using? If GDAL can't read it through python, it won't read it neither with the command line utilities.
      Maybe you can transform your file into GeoTIFF or any other available format. Usually, ENVI format is an easy one.

      Delete
  4. Hello Roger,

    I am using a SAR CEOS data file which comes with leader file, volume directory file and null volume directory file.

    I am opening a data file using gdal ReadasArray and doing operations on that 2d array,after that I want to save this 2d Array as an ENVI binary file.

    Kindly guide how to do this.

    ReplyDelete
  5. Nice blog, keep more updates about this type of information. Visit for the best Website Designing and Development Company in Delhi.
    SEO Service in Delhi

    ReplyDelete
  6. Get Mutual Fund Investment Schemes by Mutual Fund Wala and know about the best investment platform for you, to get profit.
    Best Performing Mutual Fund

    ReplyDelete
  7. http://www.duasinislam.com/pari-ka-amal/pari-ko-pane-ka-amal/

    ReplyDelete
  8. This is really nice, thanks again for this wonderful and valuable information sharing with us. Visit Kalakutir for Godown Floor Marking Painting and Base Company Logo Painting.
    Base Company Logo Painting

    ReplyDelete
  9. Thenutritionsclinic Even though we're surrounded by way of hyper speedy sports, we grow to be uninterested inside the possibilities surrounding us. This definitely is not a length to be ignored, as in this emotional country we miss all of the fortunes that come our manner. This loss of hobby might also come in our early young adults or in mid age, which we frequently call as mid life crisis. We find that we've abruptly lost hobby inside the happenings around us and hold the whole thing to ourselves.
    https://thenutritionsclinic.com/

    ReplyDelete
  10. ketogenicpedia Check out these brief hints that might even keep your existence. Short hints which can beautify your health and save your life wash fruits and vegetables earlier than you prepare dinner them. Regardless of how "natural" you located your food is, you may never actually make sure besides you in my view increase it.
    https://ketogenicpedia.com/

    ReplyDelete
  11. Nice explanation specially of the raster calculation. But I really need someone who will guide me through the path.

    Pop singer

    ReplyDelete
  12. I just loved your article on the beginners guide to starting a blog.If somebody take this blog article seriously in their life, he/she can earn his living by doing blogging.thank you for thizs article. python online training

    ReplyDelete
  13. Thanks for these great instructions. I wondered whether there is another way to import the functions since Spyder shows me warnings all the time because it doesn't know the functions imported with the asterisks were defined. Is there a way to import them directly (I couldn't find from where they are)?

    ReplyDelete
  14. To increase your visitor towards your business, you must have the best website designing company in Gurgaon which can provide you the best website design for your offline business. Website is the first impression of your business nowadays.

    ReplyDelete
  15. Interested to buy an iPhone Insurance? If you are an owner of a brand new iPhone or want to be owner then you must consider the best iPhone insurance. Connect with Your Mobile Insurance to get the best insurance for your phone.

    ReplyDelete
  16. for best quotes images in hindi to send to your family members, friends, family, gf/bf you can check out : Best HD quotes images in hindi



    ReplyDelete
  17. Are you looking for the best Website design company in Dubai? If yes then connect with CSS Founder that provides always best website design & development services at an affordable price.

    ReplyDelete
  18. Great post you shared. Thank you for this and keep posting & sharing like this. Residential projects in sector 79 Noida

    ReplyDelete
  19. Looking for Sun Shades that can help you. Connect with Al Mumtaz Tents, we can be the right choice for you.

    ReplyDelete
  20. Do you want to grow your online business? if yes then you must hire a best SEO Company in thane that can help you in boosting your business as a brand.

    ReplyDelete
  21. Have you understood the importance of websites in your business? When it comes to online business, you want to avoid shortcuts. We have put the results of a beautifully designed robust website, would you like to expand your business through the website, contact us.
    Web design company Atlantis

    ReplyDelete
  22. I got something valuable information from your blog. Thanks for posting it!
    Vastu Consultant for Business and Growth in Dubai

    ReplyDelete
  23. I would like to say that this blog really convinced me to do it! Thanks, very good post. Website Design Company in moga

    ReplyDelete
  24. Interesting post. I Have Been wondering about this issue, so thanks for posting. Pretty cool post.It 's really very nice and Useful post.Thanks Resorts in Lansdowne

    ReplyDelete
  25. If you want to get SEO of website or website, then contact with CSS founder only. Because of website design and website SEO, no one can do this work better than CSS Founder.
    Website design company sofia

    ReplyDelete
  26. With the help of Google ADS, you can give a new identity to your business and you can easily reach your business many people, for this, you have to contact Promotedial. Google Ads Services in Siliguri

    ReplyDelete
  27. 우리카지노는 대한민국의 바카라 업계를 장악하고 있는 카지노사이트 입니다. 우리카지노가 대한 민국에서 장악한 바카라 시장점유율이 50%가 넘고 10 년 넘게 온라인 바카라 시장을 장악해왔기
    때문에 대한민국에서는 우리카지노를 모르는 사람은 드뭅니다. 이런 바카라 업계의 독보적인 입지 때문에 늘 유명하거나 최고만을 찾는 사람들이 카지노사이트를 찾을때는 늘 우리카지노를 찾습니다.바카라를 처음 시작하시는 초보자분들에게도 우리카지노에서 카지노사이트를 시작하시기 좋은 환경입니다. 우리카지노사이트에서는 신규가입시 3 만쿠폰을 지급 해주기 때문입니다. 사람들이 늘 1 등만을 찾는 이유는 분명 있습니다. 다른 카지노사이트와는 달리 우리카지노를 이용하실시 에이전트를 끼고 게임을 하신다면 본사 이외에 활동쿠폰 및 오링쿠폰을 별도로 제공해주고 있기 때문입니다. 이러한 이유들 때문에 카지노사이트
    업계에서 바카라를 즐기신다면 다들 우리카지노를 선호 하십니다.
    카지노사이트에서 바카라를 이기기 물론 어렵습니다. 하지만 우리카지노의 에이전트를 끼고 바카라를 즐기신다면 승산이 있다고 봅니다. 우리카지노 에이전트의 연락처는 홈페이지로 연락하시면 언제든지 부담없이 소통가능 합니다. 카지노사이트를 선정할때는 바카라를 다른곳보다 유리하게 즐길 수 있는 카지노를 선택해야한다고 생각합니다. 그것이 바로 우리카지노 입니다. 이상으로 우리카지노와 바카라 카지노사이트 사이의 상관관계를 알아보았습니다바카라사이트.

    ReplyDelete
  28. Get a high ranking with Promotedial we are the best SEO company in SEO company in Doha . If you are looking for an SEO service agency so connect with us. 

    ReplyDelete
  29. Hurry up for the best offer for Python Training in Chennai And also get a great idea with other courses like Cyber Security, Graphic Design and Animation, Block Security, Java, Cyber Security, Oracle, Big Data, Azure, Manual and Automation Testing, DevOps, Medical Coding etc., and we also provide best technical trainers with excellent training 100+ Live Practical Sessions with Real-Time scenarios at the end of the course the freshers, experienced, and Tech professionals will be able to obtain more knowledge of the course and be able to get through the interviews on top MNC’s with an amazing package. For more details ring us up on 7504633633, 7502633633.

    ReplyDelete
  30. Thank you because you have been willing to share information with us. we will always appreciate all you have done here Website design company in Bhopal

    ReplyDelete
  31. We are the best SEO company in Giza Egypt If you are looking for an SEO company to increase your website ranking so connect with us. We can provide you SEO service at an affordable Price.

    ReplyDelete
  32. Thanks a lot . "Page load speed is an important Google ranking factor.” Seo company in Perth

    ReplyDelete
  33. “There are three responses to a piece of design – yes, no, and WOW! Wow is the one to aim for” website design company in essen

    ReplyDelete
  34. “We wander for distraction, but we travel for fulfilment.” Turkey visa

    ReplyDelete
  35. Not too lengthy. Crisp and the meaning is conveyed. Nice blog.

    Generator on rent

    ReplyDelete
  36. Heroes Blogging: Government scheme, Business Ideas, Bollywood News, Computer, Marketing Tips, Travel how to earn money and useful information!

    ReplyDelete
  37. Get your visa at an affordable price. You can easily avail it on India Visa. We are a government-certified visa company and ensure to provide your visa at a given period of time.

    ReplyDelete
  38. CSS Founder is a well-known and reputed website designing company in Delhi. With a dedicated team of experts, they provide top-notch and innovative website design services to their clients. Whether it's a small business or a large enterprise, CSS Founder has the expertise to meet the unique needs of each client, resulting in visually stunning and user-friendly websites. Their commitment to delivering quality and timely solutions has earned them the trust and loyalty of their customers. CSS Founder is truly a reliable and trustworthy option for businesses that want to establish a strong online presence and increase their digital success.

    ReplyDelete