diff --git a/plio/spatial/footprint.py b/plio/spatial/footprint.py
new file mode 100644
index 0000000000000000000000000000000000000000..43ac6490ee4f23bd2358f9928690868345c25456
--- /dev/null
+++ b/plio/spatial/footprint.py
@@ -0,0 +1,59 @@
+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
diff --git a/plio/utils/generate_vrt.py b/plio/utils/generate_vrt.py
new file mode 100644
index 0000000000000000000000000000000000000000..e58d755d0767cf28fc2817fe373f6a5613204065
--- /dev/null
+++ b/plio/utils/generate_vrt.py
@@ -0,0 +1,48 @@
+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)