Merge pull request #65 from jlaura/camera

Adds CSM Camera Module

plio/camera/csm.py

+import datetime
+import json
+
+try:
+    import usgscam as cam
+    from cycsm.isd import Isd
+    camera_support = True
+except:
+    camera_support = False
+import requests
+
+from plio.utils.utils import find_in_dict
+from plio.io.io_json import NumpyEncoder
+
+
+def data_from_cube(header):
+    data = {}
+    data['START_TIME'] = find_in_dict(header, 'StartTime')
+    data['SPACECRAFT_NAME'] = find_in_dict(header, 'SpacecraftName')
+    data['INSTRUMENT_NAME'] = find_in_dict(header, 'InstrumentId')
+    data['SAMPLING_FACTOR'] = find_in_dict(header, 'SpatialSumming')
+    data['SAMPLE_FIRST_PIXEL'] = find_in_dict(header, 'SampleFirstPixel')
+    data['IMAGE'] = {}
+    data['IMAGE']['LINES'] = find_in_dict(header, 'Lines')
+    data['IMAGE']['SAMPLES'] = find_in_dict(header, 'Samples')
+    data['TARGET_NAME'] = find_in_dict(header, 'TargetName')
+    data['LINE_EXPOSURE_DURATION'] = find_in_dict(header, 'LineExposureDuration')
+    data['SPACECRAFT_CLOCK_START_COUNT'] = find_in_dict(header, 'SpacecraftClockCount')
+    return data
+
+def create_camera(obj, url='http://smalls:8002/api/1.0/missions/mars_reconnaissance_orbiter/csm_isd'):
+    if not camera_support:
+        print("Usgscam library not installed. Camera capabilities are disabled")
+    
+    data = json.dumps(data_from_cube(obj.metadata), cls=NumpyEncoder)
+    r = requests.post(url, data=data)
+
+    # Get the ISD back and instantiate a local ISD for the image
+    isd = r.json()['data']['isd']
+    i = Isd.loads(isd)
+
+    # Create the plugin and camera as usual
+    plugin = cam.genericls.Plugin()
+    return plugin.from_isd(i, plugin.modelname(0))
\ No newline at end of file

plio/camera/tests/test_csm.py

+import datetime
+import json
+from unittest import mock
+
+import pytest
+import pvl
+
+import usgscam
+
+from plio.camera import csm
+from plio.examples import get_path
+
+def mock_requests_post(*args, **kwargs):
+    class MockResponse:
+        def __init__(self, json_data, status_code):
+            self.json_data = json_data
+            self.status_code = status_code
+
+        def json(self):
+            return self.json_data
+
+    if 'mars_reconnaissance_orbiter' in args[0]:
+        with open(get_path('ctx.response'), 'r') as f:
+            resp = json.load(f)
+        return MockResponse(resp, 200)
+
+    return MockResponse(None, 404)
+
+
+@pytest.fixture
+def header():
+    return pvl.load(get_path('ctx.pvl'))
+
+@pytest.fixture
+def req_obj():
+    return 
+
+def test_data_from_cube(header):
+    data = csm.data_from_cube(header)
+    assert data['START_TIME'] == datetime.datetime(2008, 9, 17, 5, 8, 10, 820000)
+
+@mock.patch('requests.post', side_effect=mock_requests_post)
+def test_create_camera(header):
+    cam = csm.create_camera(header)
+    assert isinstance(cam, usgscam.genericls.SensorModel)

plio/io/io_json.py

 import json
+import numpy as np
 
 
+class NumpyEncoder(json.JSONEncoder):
+    def default(self, obj):
+        if isinstance(obj, np.ndarray):
+            return obj.tolist()
+        elif isinstance(obj, datetime.date):
+            return obj.isoformat()
+        return json.JSONEncoder.default(self, obj)
+
 def read_json(inputfile):
     """
     Read the input json file into a python dictionary.

plio/spatial/footprint.py

+import numpy as np
+import pyproj
+import ogr
+
+def generate_gcps(camera, nnodes=5, semi_major=3396190, semi_minor=3376200):
+    ecef = pyproj.Proj(proj='geocent', a=semi_major, b=semi_minor)
+    lla = pyproj.Proj(proj='latlon', a=semi_major, b=semi_minor)
+
+    isize = camera.imagesize[::-1]
+    x = np.linspace(0,isize[1], 10)
+    y = np.linspace(0,isize[0], 10)
+    boundary = [(i,0.) for i in x] + [(isize[1], i) for i in y[1:]] +\
+               [(i, isize[0]) for i in x[::-1][1:]] + [(0.,i) for i in y[::-1][1:]]
+    gnds = np.empty((len(boundary), 3))
+    for i, b in enumerate(boundary):
+        gnds[i] = camera.imageToGround(*b, 0)
+    lons, lats, alts = pyproj.transform(ecef, lla, gnds[:,0], gnds[:,1], gnds[:,2])
+    lla = np.vstack((lons, lats, alts)).T
+
+    tr = zip(boundary, lla)
+
+    gcps = []
+    for i, t in enumerate(tr):
+        l = '<GCP Id="{}" Info="{}" Pixel="{}" Line="{}" X="{}" Y="{}" Z="{}" />'.format(i, i, t[0][1], t[0][0], t[1][0], t[1][1], t[1][2])
+        gcps.append(l)
+
+    return gcps
+
+def generate_latlon_footprint(camera, nnodes=5, semi_major=3396190, semi_minor=3376200):
+    ecef = pyproj.Proj(proj='geocent', a=semi_major, b=semi_minor)
+    lla = pyproj.Proj(proj='latlon', a=semi_major, b=semi_minor)
+
+    isize = camera.imagesize[::-1]
+    x = np.linspace(0,isize[1], 10)
+    y = np.linspace(0,isize[0], 10)
+    boundary = [(i,0.) for i in x] + [(isize[1], i) for i in y[1:]] +\
+               [(i, isize[0]) for i in x[::-1][1:]] + [(0.,i) for i in y[::-1][1:]]
+    ring = ogr.Geometry(ogr.wkbLinearRing)
+    for i in boundary:
+        gnd = camera.imageToGround(*i, 0)
+        lons, lats, alts = pyproj.transform(ecef, lla, gnd[0], gnd[1], gnd[2])
+        ring.AddPoint(lons, lats)
+    poly = ogr.Geometry(ogr.wkbPolygon)
+    poly.AddGeometry(ring)
+    return poly
+
+def generate_bodyfixed_footprint(camera, nnodes=5):
+    isize = camera.imagesize[::-1]
+    x = np.linspace(0,isize[1], 10)
+    y = np.linspace(0,isize[0], 10)
+    boundary = [(i,0.) for i in x] + [(isize[1], i) for i in y[1:]] +\
+               [(i, isize[0]) for i in x[::-1][1:]] + [(0.,i) for i in y[::-1][1:]]
+    ring = ogr.Geometry(ogr.wkbLinearRing)
+    for i in boundary:
+        gnd = camera.imageToGround(*i, 0)
+        ring.AddPoint(gnd[0], gnd[1], gnd[2])
+    poly = ogr.Geometry(ogr.wkbPolygon)
+    poly.AddGeometry(ring)
+    return poly

plio/utils/generate_vrt.py

+import gdal
+import os
+import jinja2
+import numpy as np
+
+from plio.camera.footprint import generate_gcps
+
+def warped_vrt(camera, raster_size, fpath, outpath=None, no_data_value=0):
+    gcps = generate_gcps(camera)
+    xsize, ysize = raster_size
+
+    if outpath is None:
+        outpath = os.path.dirname(fpath)
+    outname = os.path.splitext(os.path.basename(fpath))[0] + '.vrt'
+    outname = os.path.join(outpath, outname)
+
+    xsize, ysize = raster_size
+    vrt = r'''<VRTDataset rasterXSize="{{ xsize }}" rasterYSize="{{ ysize }}">
+     <Metadata/>
+     <GCPList Projection="{{ proj }}">
+     {% for gcp in gcps -%}
+       {{gcp}}
+     {% endfor -%}
+    </GCPList>
+     <VRTRasterBand dataType="Float32" band="1">
+       <NoDataValue>{{ no_data_value }}</NoDataValue>
+       <Metadata/>
+       <ColorInterp>Gray</ColorInterp>
+       <SimpleSource>
+         <SourceFilename relativeToVRT="0">{{ fpath }}</SourceFilename>
+         <SourceBand>1</SourceBand>
+         <SourceProperties rasterXSize="{{ xsize }}" rasterYSize="{{ ysize }}"
+    DataType="Float32" BlockXSize="512" BlockYSize="512"/>
+         <SrcRect xOff="0" yOff="0" xSize="{{ xsize }}" ySize="{{ ysize }}"/>
+         <DstRect xOff="0" yOff="0" xSize="{{ xsize }}" ySize="{{ ysize }}"/>
+       </SimpleSource>
+     </VRTRasterBand>
+    </VRTDataset>'''
+
+    context = {'xsize':xsize, 'ysize':ysize,
+               'gcps':gcps,
+               'proj':'+proj=longlat +a=3396190 +b=3376200 +no_defs',
+               'fpath':fpath,
+               'no_data_value':no_data_value}
+    template = jinja2.Template(vrt)
+    tmp = template.render(context)
+    warp_options = gdal.WarpOptions(format='VRT', dstNodata=0)
+    gdal.Warp(outname, tmp, options=warp_options)

recipe/meta.yaml

@@ -42,6 +42,8 @@ requirements:
     - scipy
     - certifi
     - pyproj
+    - usgscam
+    - jinja2
 
 test:
   imports: