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
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)
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
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
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)
instancing a satellite data object, getting numpy arrays from single channels (a part of the example/products mentioned early)
featuring python, numpy, gdal-python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 | #!/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... |