Merge d1bf1b6 into 36c3e4e

etc/s2a.cfg.template

@@ -0,0 +1,88 @@
+[satellite]
+satname = s2a
+variant = 
+number = 
+instruments = ('msi',)
+
+[msi-level2]
+filename=%(satellite)s_OPER_%(instrument)s_L1C_??_????_%Y%m%dT%H%M%S_A%(orbit)s_T?????_%(band)s.jp2
+dir = /s2a_datadir
+format = s2_msi
+
+[msi-b1]
+frequency = (0.433, 0.443, 0.453)
+resolution = 60
+name = 'B01'
+size = (1830,1830)
+
+[msi-b2]
+frequency = (0.4575, 0.490, 0.5225)
+resolution = 10
+name = 'B02'
+size = (10980,10980)
+
+[msi-b3]
+frequency = (0.5225, 0.560, 0.5775)
+resolution = 10
+name = 'B03'
+size = (10980,10980)
+
+[msi-b4]
+frequency = (0.650, 0.665, 0.680)
+resolution = 10
+name = 'B04'
+size = (10980,10980)
+
+[msi-b5]
+frequency = (0.6975, 0.705, 0.7125)
+resolution = 20
+name = 'B05'
+size = (5490,5490)
+
+[msi-b6]
+frequency = (0.7325, 0.740, 0.7475)
+resolution = 20
+name = 'B06'
+size = (5490,5490)
+
+[msi-b7]
+frequency = (0.773, 0.783, 0.793)
+resolution = 20 
+name = 'B07'
+size = (5490,5490)
+
+[msi-b8]
+frequency = (0.7845, 0.842, 0.8995)
+resolution = 10  
+name = 'B08'
+size = (10980,10980)
+
+[msi-b8a]
+frequency = (0.855,0.865,0.875)
+resolution = 20
+name = 'B8A'
+size = (5490,5490)
+
+[msi-b9]
+frequency = (0.935,0.945,0.955)
+resolution = 60   
+name = 'B09'
+size = (1830,1830)
+
+[msi-b10]
+frequency = (1.360,1.375,1.390)
+resolution = 60
+name = 'B10'
+size = (1830,1830)
+
+[msi-b11]
+frequency = (1.565,1.610,1.655)
+resolution = 20
+name = 'B11'
+size = (5490,5490)
+
+[msi-b12]
+frequency = (2.100,2.190,2.280)
+resolution = 20 
+name = 'B12'
+size = (5490,5490)

mpop/instruments/s2_composites.py

@@ -0,0 +1,24 @@
+from mpop.imageo.geo_image import GeoImage
+
+def s2_truecolor(self):
+
+	self.check_channels('B02','B03','B04')
+
+	ch1 = self['B04'].data
+	ch2 = self['B03'].data
+	ch3 = self['B02'].data
+
+	img = GeoImage((ch1, ch2, ch3),
+				self.area,
+                                 self.time_slot,
+                                 fill_value=None,
+                                 mode="RGB")
+
+	img.enhance(stretch="linear")
+	#img.enhance(stretch="histogram")
+	img.enhance(gamma=2.0)
+
+	return img
+
+s2_truecolor.prerequisites = set(['B02', 'B03','B04'])
+msi=[s2_truecolor]

mpop/satin/s2_msi.py

@@ -0,0 +1,60 @@
+#!/usr/bin/python
+"""Loader for s2, jpeg2000 format.
+"""
+#Matias Takala FMI 2016
+
+from ConfigParser import ConfigParser
+from mpop import CONFIG_PATH
+import glymur
+import os
+from pyresample import utils
+import pickle,re
+from mpop.satellites import GenericFactory
+import glob
+import numpy.ma as ma
+
+#in this version Q_V is hardcoded but could be read from metadata
+QUANTIFICATION_VALUE=10000
+
+def parse_tile(file):
+	tile=re.findall('T(\d{2}\w{3})_',file)	
+	f=open('s2tiles.pickle','r')
+	s2tiles=pickle.load(f)
+	f.close()
+	return [tile[0],s2tiles[tile[0]]]
+
+def read_jp2_data(file):
+	jp2=glymur.Jp2k(file)
+	data=jp2[:]/(QUANTIFICATION_VALUE+0.0)
+	return data
+
+def open_s2_tile(fname):
+	data=read_jp2_data(fname)	
+	size=data.shape
+	params=parse_tile(fname)
+	areadef=utils.get_area_def(params[0],"Sentinel 2 tile "+params[0],'PROJ EPSG:'+params[1][0],'init=epsg:'+params[1][0],size[0],size[1],params[1][1])
+	return([data,areadef])
+
+def load(satscene):
+    """Load jpeg2000 data.
+    """
+
+    # Read config file content
+    conf = ConfigParser()
+    conf.read(os.path.join(CONFIG_PATH, satscene.fullname + ".cfg"))
+
+    for chn_name in satscene.channels_to_load:
+	    values = {"orbit": satscene.orbit,
+       		      "satname": satscene.satname.upper(),
+	              "number": satscene.number,
+	              "instrument": satscene.instrument_name.upper(),
+	              "satellite": satscene.fullname.upper(),
+		      "band": chn_name		     	
+	              }
+	    filename = os.path.join(
+	    conf.get("msi-level2", "dir"),
+       	    satscene.time_slot.strftime(conf.get("msi-level2","filename",raw=True)) % values) 
+	    filelist=glob.glob(filename)
+	    data_area=open_s2_tile(filelist[0])
+	    satscene[chn_name] = ma.masked_array(data_area[0])
+	    satscene[chn_name].area=data_area[1]

utils/get_tile_def.py

@@ -0,0 +1,51 @@
+#!/usr/bin/python
+
+import xml.etree.ElementTree as ET
+from pyresample import utils
+import pickle
+import urllib2
+
+length=109800
+
+#https://sentinel.esa.int/documents/247904/1955685/S2A_OPER_GIP_TILPAR_20150622T000000_21000101T000000_ZZ_0001
+FNAME="S2A_OPER_GIP_TILPAR_20150622T000000_21000101T000000_ZZ_0001"
+
+TILPAR_URL="https://sentinel.esa.int/documents/247904/1955685/"+FNAME
+
+FNAME=FNAME+".kml"
+
+tiles = urllib2.urlopen(TILPAR_URL)
+with open(FNAME,'wb') as output:
+  output.write(tiles.read())
+
+tiles.close()
+
+tree = ET.parse(FNAME)
+root = tree.getroot()
+
+s2tiles={}
+
+for pm in root.iter('{http://www.opengis.net/kml/2.2}Placemark'):
+	tilename=None
+	epsg=None
+	utm_ul_x=None
+	utm_ul_y=None
+
+	for name in pm.iter('{http://www.opengis.net/kml/2.2}name'):
+		tilename=name.text
+	for simple in pm.iter('{http://www.opengis.net/kml/2.2}SimpleData'):
+		if (simple.attrib['name']=='epsg'):
+			epsg=simple.text
+		if(simple.attrib['name']=='utm_ul_x'):
+			utm_ul_x=simple.text
+		if(simple.attrib['name']=='utm_ul_y'):
+			utm_ul_y=simple.text
+
+	extent=(float(utm_ul_x),float(utm_ul_y)-length,float(utm_ul_x)+length,float(utm_ul_y))
+
+	s2tiles[tilename]=[epsg,extent]
+
+f=open('s2tiles.pickle','w')
+pickle.dump(s2tiles,f)
+f.close()
+