Sean Gillies (Posts about warp)https://sgillies.net/tags/warp.atom2023-12-31T01:26:24ZSean GilliesNikolaWarping images with rasteriohttps://sgillies.net/2014/02/25/warping-images-with-rasterio.html2014-02-25T00:00:00-07:002014-02-25T00:00:00-07:00Sean Gillies<p>Warping imagery from one spatial basis to another is a crucial function. I've
been determined that rasterio be able to do this. Scipy provides
a <a class="reference external" href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.interpolation.geometric_transform.html#scipy.ndimage.interpolation.geometric_transform">geometric_transform()</a>
function that I like, but I'd like rasterio to operate at a higher level and
take advantage of GDAL's fast transformers. GDAL's gdalwarp program only
operates on files, and I'm aiming lower than this (so rasterio can let us build
better scripts). The <code class="docutils literal">ReprojectImage()</code> function in osgeo.gdal, as pointed
out in <a class="reference external" href="http://jgomezdans.github.io/gdal_notes/reprojection.html">http://jgomezdans.github.io/gdal_notes/reprojection.html</a>, only works on
GDAL dataset objects. And the code to make them gets messy.</p>
<p>My requirements for a rasterio warper are:</p>
<ul class="simple">
<li><p>Speed approaching gdalwarp.</p></li>
<li><p>Same results (same algorithms) as gdalwarp.</p></li>
<li><p>Sources and destinations can be ndarrays as well as bands in files.</p></li>
<li><p>Simplicity and clarity.</p></li>
</ul>
<p>I think I've managed to satisfy these requirements. The warp module will be the
key feature of the next rasterio release.</p>
<p>Here's an example within the interactive dataset inspector.</p>
<div class="code"><pre class="code console"><a id="rest_code_422ba25feb1343ed89ffc684f1055c17-1" name="rest_code_422ba25feb1343ed89ffc684f1055c17-1" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_422ba25feb1343ed89ffc684f1055c17-1"></a><span class="gp">$ </span>rasterio.insp rasterio/tests/data/RGB.byte.tif
<a id="rest_code_422ba25feb1343ed89ffc684f1055c17-2" name="rest_code_422ba25feb1343ed89ffc684f1055c17-2" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_422ba25feb1343ed89ffc684f1055c17-2"></a><span class="go">Rasterio 0.7 Interactive Inspector (Python 2.7.5)</span>
<a id="rest_code_422ba25feb1343ed89ffc684f1055c17-3" name="rest_code_422ba25feb1343ed89ffc684f1055c17-3" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_422ba25feb1343ed89ffc684f1055c17-3"></a><span class="go">Type "src.name", "src.read_band(1)", or "help(src)" for more information.</span>
</pre></div>
<p>The <code class="docutils literal">reproject()</code> function takes source and destination ndarrays and an
affine transformation vector and coordinate reference system for each. In
this example, I'm warping the UTM source to a Web Mercator destination of slightly
smaller extent. I computed the extents of the destination using pyplot, rasterio
has yet no method to come up with an optimal destination extent (to do).</p>
<div class="code"><pre class="code pycon"><a id="rest_code_f1a546b265df48539af28051a906edcf-1" name="rest_code_f1a546b265df48539af28051a906edcf-1" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-1"></a><span class="gp">>>> </span><span class="n">source</span> <span class="o">=</span> <span class="n">src</span><span class="o">.</span><span class="n">read_band</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-2" name="rest_code_f1a546b265df48539af28051a906edcf-2" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-2"></a><span class="gp">>>> </span><span class="n">dst_transform</span> <span class="o">=</span> <span class="p">[</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-3" name="rest_code_f1a546b265df48539af28051a906edcf-3" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-3"></a><span class="gp">... </span> <span class="o">-</span><span class="mf">8789636.708</span><span class="p">,</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-4" name="rest_code_f1a546b265df48539af28051a906edcf-4" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-4"></a><span class="gp">... </span> <span class="mf">300.0</span><span class="p">,</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-5" name="rest_code_f1a546b265df48539af28051a906edcf-5" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-5"></a><span class="gp">... </span> <span class="mf">0.0</span><span class="p">,</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-6" name="rest_code_f1a546b265df48539af28051a906edcf-6" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-6"></a><span class="gp">... </span> <span class="mf">2943560.235</span><span class="p">,</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-7" name="rest_code_f1a546b265df48539af28051a906edcf-7" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-7"></a><span class="gp">... </span> <span class="mf">0.0</span><span class="p">,</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-8" name="rest_code_f1a546b265df48539af28051a906edcf-8" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-8"></a><span class="gp">... </span> <span class="o">-</span><span class="mf">300.0</span><span class="p">]</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-9" name="rest_code_f1a546b265df48539af28051a906edcf-9" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-9"></a><span class="gp">>>> </span><span class="n">dst_crs</span> <span class="o">=</span> <span class="p">{</span><span class="s1">'init'</span><span class="p">:</span> <span class="s1">'EPSG:3857'</span><span class="p">}</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-10" name="rest_code_f1a546b265df48539af28051a906edcf-10" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-10"></a><span class="gp">>>> </span><span class="n">destin</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">empty</span><span class="p">(</span><span class="n">src</span><span class="o">.</span><span class="n">shape</span><span class="p">,</span> <span class="n">dtype</span><span class="o">=</span><span class="n">np</span><span class="o">.</span><span class="n">uint8</span><span class="p">)</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-11" name="rest_code_f1a546b265df48539af28051a906edcf-11" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-11"></a><span class="gp">>>> </span><span class="kn">from</span> <span class="nn">rasterio</span> <span class="kn">import</span> <span class="n">_warp</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-12" name="rest_code_f1a546b265df48539af28051a906edcf-12" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-12"></a><span class="gp">>>> </span><span class="n">_warp</span><span class="o">.</span><span class="n">reproject</span><span class="p">(</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-13" name="rest_code_f1a546b265df48539af28051a906edcf-13" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-13"></a><span class="gp">... </span> <span class="n">source</span><span class="p">,</span> <span class="n">destin</span><span class="p">,</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-14" name="rest_code_f1a546b265df48539af28051a906edcf-14" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-14"></a><span class="gp">... </span> <span class="n">src</span><span class="o">.</span><span class="n">transform</span><span class="p">,</span> <span class="n">src</span><span class="o">.</span><span class="n">crs</span><span class="p">,</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-15" name="rest_code_f1a546b265df48539af28051a906edcf-15" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-15"></a><span class="gp">... </span> <span class="n">dst_transform</span><span class="p">,</span> <span class="n">dst_crs</span><span class="p">)</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-16" name="rest_code_f1a546b265df48539af28051a906edcf-16" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-16"></a><span class="gp">>>> </span><span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-17" name="rest_code_f1a546b265df48539af28051a906edcf-17" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-17"></a><span class="gp">>>> </span><span class="n">plt</span><span class="o">.</span><span class="n">subplot</span><span class="p">(</span><span class="mi">121</span><span class="p">)</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-18" name="rest_code_f1a546b265df48539af28051a906edcf-18" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-18"></a><span class="go"><matplotlib.axes.AxesSubplot object at 0x1075362d0></span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-19" name="rest_code_f1a546b265df48539af28051a906edcf-19" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-19"></a><span class="gp">>>> </span><span class="n">plt</span><span class="o">.</span><span class="n">imshow</span><span class="p">(</span><span class="n">source</span><span class="p">)</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-20" name="rest_code_f1a546b265df48539af28051a906edcf-20" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-20"></a><span class="go"><matplotlib.image.AxesImage object at 0x1078b1450></span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-21" name="rest_code_f1a546b265df48539af28051a906edcf-21" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-21"></a><span class="gp">>>> </span><span class="n">plt</span><span class="o">.</span><span class="n">gray</span><span class="p">()</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-22" name="rest_code_f1a546b265df48539af28051a906edcf-22" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-22"></a><span class="gp">>>> </span><span class="n">plt</span><span class="o">.</span><span class="n">subplot</span><span class="p">(</span><span class="mi">122</span><span class="p">)</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-23" name="rest_code_f1a546b265df48539af28051a906edcf-23" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-23"></a><span class="go"><matplotlib.axes.AxesSubplot object at 0x1075ffe50></span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-24" name="rest_code_f1a546b265df48539af28051a906edcf-24" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-24"></a><span class="gp">>>> </span><span class="n">plt</span><span class="o">.</span><span class="n">imshow</span><span class="p">(</span><span class="n">destin</span><span class="p">)</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-25" name="rest_code_f1a546b265df48539af28051a906edcf-25" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-25"></a><span class="go"><matplotlib.image.AxesImage object at 0x1078b19d0></span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-26" name="rest_code_f1a546b265df48539af28051a906edcf-26" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-26"></a><span class="gp">>>> </span><span class="n">plt</span><span class="o">.</span><span class="n">gray</span><span class="p">()</span>
<a id="rest_code_f1a546b265df48539af28051a906edcf-27" name="rest_code_f1a546b265df48539af28051a906edcf-27" href="https://sgillies.net/2014/02/25/warping-images-with-rasterio.html#rest_code_f1a546b265df48539af28051a906edcf-27"></a><span class="gp">>>> </span><span class="n">plt</span><span class="o">.</span><span class="n">savefig</span><span class="p">(</span><span class="s1">'/tmp/warp.png'</span><span class="p">)</span>
</pre></div>
<p>It's pretty much just three statements. One to get a source ndarray from a
dataset. One to create a destination ndarray. One to call <code class="docutils literal">reproject()</code>.</p>
<p>I'm out of time to properly transform the axes in the plots to show coordinates
in their proper CRS, but they're warped correctly.</p>
<img alt="http://farm4.staticflickr.com/3776/12773733534_c51761076e_z_d.jpg" src="http://farm4.staticflickr.com/3776/12773733534_c51761076e_z_d.jpg">