Examples - ARPA-SIMC/meteosatlib GitHub Wiki

accessing datasets

when using meteosatlib command line tools it's possible to access a dataset both by pointing at one of its files:

msat --view /somepath/H-000-MSG3__-MSG3________-IR_039___-000005___-201307100230-C_

or by using the special syntax [RESOLUTION]:[SATELLITE]:[CHANNEL]:[TIMESTAMP]

msat --view /somepath/H:MSG3:IR_039:201307100230

the latter method allows to access the reflectance and solar zenith angle special datasets:

msat --view /somepath/H:MSG3:IR_039r:201307100230

msat --view /somepath/H:MSG3:IR_039a:201307100230

Note:

• the reflectance dataset is available since meteosatlib 1.2 and only for VIS006 VIS008 IR_016 IR_039 channels
• the solar zenith angle dataset is available since meteosatlib 1.4*

cropping

graphical tool for finding out boundary box vertices:

msat-view H:MSG2:HRV:201001191215

then:

# area defined in pixel
msat   --display  --area="5760,1792,1126,1254" H:MSG2:HRV:201203130715
# lat/lon defined area, the same of the above example
msat   --display  --Area="31.5872,48.1860,2.7367,22.2706" H:MSG2:HRV:201203130715

(see manpage for details on area cropping options)

output formats

as the manpage says:

   -c, --conv FMT
          Convert to the given format (see gdalinfo --formats for a list)

   --jpg  Convert to JPEG (with gray scale normalization).

   --png  Convert to PNG (with gray scale normalization).

meaning you can have two (very) different kinds of graphic output

# gray scale normalized png (increased contrast):
msat --png --area="5760,1792,1126,1254" H:MSG2:HRV:201203130715
# png:
msat -c PNG --area="5760,1792,1126,1254" H:MSG2:HRV:201203130715

plus anything that gdal can handle

# grib:
msat -c MsatGRIB   --area="5056,2000,448,2000" H:MSG2:WV_073:201001191215
# netcdf:
msat -c MsatNetCDF --area="5056,2000,448,2000" H:MSG2:WV_073:201001191215

reprojecting data

this is more like a gdal feature, but you can use it on meteosat data thanks to msat-gdal driver, so here it is:

# to lat/lon geotiff
gdalwarp -t_srs "+proj=latlong" -te -10 30 30 60 -of GTiff /mydatapath/H:MSG2:HRV:201203130715 something.gtiff
# to utm32 WGS84 ENVI data format
gdalwarp -t_srs '+proj=utm +zone=32 +datum=WGS84' -te -1349315 3475750 1657788 6838296 -of ENVI /mydatapath/H:MSG2:IR_O16:201203130715

georeferenced extents of output file (-te option) is not mandatory but highly recommended since original meteosat data may be huge

composite RGB products

there's a python script in the meteosatlib sources (example/products) that covers part of the products described in Eumetsat's MSG Channels Interpretation Guide featuring cropping, reprojection, shapefile overlaying and a help, too:

$ ./products --help
Usage: products [options]

Generate satellite products

Options:
  --version             show program's version number and exit
  -h, --help            show this help message and exit
  -q, --quiet           quiet mode: only output fatal errors
  -v, --verbose         verbose mode
  -s SRCDIR, --srcdir=SRCDIR
                        directory with the HRIT data. Default:
                        /autofs/scratch1/satope/done/
  -t TIME, --time=TIME  datetime, as 'YYYYMMDDHHMM', default: 201212121500
  -a x,dx,y,dy, --area=x,dx,y,dy
                        datetime, as 'YYYYMMDDHHMM', default:
                        1350,1400,100,800
  --shp=shapefile       shapefile to use for coastlines
  -f fmt, --format=fmt  output file format. See gdalinfo --formats. Default:
                        GTiff
  -d dir, --destdir=dir
                        output directory. Default: .
  --warp=opts           gdalwarp options to use to warp input channels before
                        using them. When used, area cropping is disabled.
                        Default: none
  --cachedir=dir        cache directory for warped channels. Default: .

example:

./products -t 201212061200 -s ~/mydata/ --shp ~/myshp/worldboundaries.shp -d /tmp/products/ --warp='-t_srs "+proj=latlong" -te -10 30 30 60' --cachedir=/tmp/products/cache/

(shapefile must be in the same projection specified in --warp options)

working with data arrays

instancing a satellite data object, getting numpy arrays from single channels (a part of the example/products mentioned early)

featuring python, numpy, gdal-python

#!/usr/bin/python
import gdal
from gdalconst import *
import datetime
import os
import os.path
import numpy

SRCDIR = "/my/msg/data/path/"
AREA = (1350, 100, 1400, 800)

class Satellite(object):
    def __init__(self, dt):
        self.dt = dt

    def gdal_dataset(self, channel_name):
        fname = "H:MSG2:%s:%s" % (
                channel_name,
                self.dt.strftime("%Y%m%d%H%M")
        )
        return gdal.Open(os.path.join(SRCDIR, fname), GA_ReadOnly)

    def get_array(self, channel):
        ds = self.gdal_dataset(channel)
        rb = ds.GetRasterBand(1)
        return rb.ReadAsArray(*AREA)

sat = Satellite(datetime.datetime(2012, 11, 6, 13, 15)
ir108 = sat.get_array("IR_108")
vis06 = sat.get_array("VIS006")
# etc...