cookie_cutter fix, add alltouch, drop py27

.circleci/config.yml

@@ -34,40 +34,6 @@ jobs:
           paths:
             - ~/project
 
-  install_27:
-    <<: *container
-    steps:  
-      - *restore_repo
-      - restore_cache:
-          keys: 
-            - v3-dependencies27-{{ checksum "requirements.txt"}}
-            - v3-dependencies27
-      - run: |
-            pip install virtualenv
-            virtualenv ~/venv27
-            . ~/venv27/bin/activate
-            pip install -r requirements.txt
-            pip install -r requirements-dev.txt
-      - save_cache:
-          key: v3-dependencies27-{{ checksum "requirements.txt"}}
-          paths:
-            - ~/venv27
-
-  test_27:
-    <<: *container
-    steps:
-      - *restore_repo
-      - restore_cache:
-          keys:
-            - v3-dependencies27-{{ checksum "requirements.txt"}}
-            - v3-dependencies27
-      - run: |
-            yum install -y swig
-            . ~/venv27/bin/activate
-            pip install -e . --no-cache
-            cd test
-            pytest --cov gippy     
-
   install_36:
     <<: *container
     steps:  
@@ -147,12 +113,3 @@ workflows:
           filters:
             branches:
               only: master
-  build_test_27:
-    jobs:
-      - checkout_code
-      - install_27:
-          requires:
-            - checkout_code
-      - test_27:
-          requires:
-            - install_27

.gitignore

@@ -14,3 +14,4 @@ gippy/gippy_wrap.cpp
 gippy/algorithms.py
 gippy/algorithms_wrap.cpp
 /.vs
+.venv*

Dockerfile

@@ -3,21 +3,19 @@ FROM developmentseed/geolambda:latest
 WORKDIR /build
 
 RUN \
-    yum install -y swig;
+    yum install -y swig \
+        && curl https://bootstrap.pypa.io/get-pip.py -o get-pip.py \
+        && python3 get-pip.py
 
 COPY requirements*txt /build/
 RUN \
-    pip2 install -r requirements.txt; \
-    pip2 install -r requirements-dev.txt; \
     pip3 install -r requirements.txt; \
     pip3 install -r requirements-dev.txt;
 
 COPY . /build
 RUN \
     git clean -xfd; \
-    pip2 install .; \
-    git clean -xfd; \
-    pip3 install .; \
+    pip3 install . ; \
     rm -rf /build/*;
 
 WORKDIR /home/geolambda

GIP/GeoAlgorithms.cpp

@@ -223,7 +223,7 @@ namespace gip {
         - if GeoFeature is provided it will it's SRS. If not, proj parameter will be used (EPSG:4326 default)
     */
     GeoImage cookie_cutter(const std::vector<GeoImage>& geoimgs, string filename,
-            GeoFeature feature, bool crop, string proj, float xres, float yres, int interpolation, dictionary options) {
+            GeoFeature feature, bool crop, string proj, float xres, float yres, int interpolation, dictionary options, bool alltouch) {
         if (Options::verbose() > 1)
             cout << "GIPPY: cookie_cutter (" << geoimgs.size() << " files) - " << filename << endl;
 
@@ -233,9 +233,12 @@ namespace gip {
             extents.push_back(i->extent());
         }
         BoundingBox ext = union_all(extents);
+        // Control special handling of extents if vector is provided
+        int vector2raster = 0 ;
 
         // if valid feature provided use that extent
         if (feature.valid()) {
+            vector2raster = 1 ;
             if (proj == "")
                 proj = feature.srs();
             // transform extent to desired srs
@@ -256,11 +259,20 @@ namespace gip {
 
         // create output
         // convert extent to resolution units
-        int xsz = std::ceil(ext.width() / std::abs(xres));
-        int ysz = std::ceil(ext.height() / std::abs(yres));
-
-        CImg<double> bbox(4,1,1,1, ext.x0(), ext.y0(), ext.width(), ext.height());
-        GeoImage imgout = GeoImage::create(filename, xsz, ysz, geoimgs[0].nbands(), 
+	// one pixel is added IFF transition from vector to raster space
+        int xsz = std::ceil(ext.width() / std::abs(xres)) + vector2raster;
+        int ysz = std::ceil(ext.height() / std::abs(yres)) + vector2raster;
+
+        double xshift = -0.5 * vector2raster * std::abs(xres);
+        double yshift = -0.5 * vector2raster * std::abs(yres);
+
+        /* Multiply the x and y size by the desired resolution to force the output
+           image to have a size evenly divisible by the res. xsz and ysz above have
+           been increased by one pixel to avoid the infamous "lost pixel"
+           when cutting a raster with a vector. Finally, to minimize pixel drift
+           amongst all of this, the whole image is shifted NW a half pixel. */
+        CImg<double> bbox(4,1,1,1, ext.x0() + xshift, ext.y0() + yshift, xsz * std::abs(xres), ysz * std::abs(yres));
+        GeoImage imgout = GeoImage::create(filename, xsz, ysz, geoimgs[0].nbands(),
                             proj, bbox, geoimgs[0].type().string(), "", false, options);
 
         imgout.add_meta(geoimgs[0].meta());
@@ -275,13 +287,13 @@ namespace gip {
         //metadata["SourceFiles"] = to_string(geoimgs.basenames());
         if (interpolation > 1) metadata["Interpolation"] = to_string(interpolation);
         imgout.add_meta(metadata);
-        
+
         bool noinit(false);
         for (unsigned int i=0; i<geoimgs.size(); i++) {
-            geoimgs[i].warp_into(imgout, feature, interpolation, noinit);
+            geoimgs[i].warp_into(imgout, feature, interpolation, noinit, alltouch);
             noinit = true;
         }
-    
+
         return imgout;
     }
 

GIP/GeoImage.cpp

@@ -326,12 +326,12 @@ namespace gip {
     }
 
 
-    GeoImage& GeoImage::warp_into(GeoImage& imgout, GeoFeature feature, int interpolation, bool noinit) const {
+    GeoImage& GeoImage::warp_into(GeoImage& imgout, GeoFeature feature, int interpolation, bool noinit, bool alltouch) const {
         if (Options::verbose() > 2) std::cout << basename() << " warping into " << imgout.basename() << std::endl;
 
         // Assume that these have the same number of bands
         for (unsigned int b=0; b<this->nbands(); b++) {
-            (*this)[b].warp_into(imgout[b], feature, interpolation, noinit);
+            (*this)[b].warp_into(imgout[b], feature, interpolation, noinit, alltouch);
         }
 
         return imgout;

GIP/GeoRaster.cpp

@@ -207,7 +207,7 @@ namespace gip {
         return *this;
     }
 
-    GeoRaster& GeoRaster::warp_into(GeoRaster& imgout, GeoFeature feature, int interpolation, bool noinit) const {
+    GeoRaster& GeoRaster::warp_into(GeoRaster& imgout, GeoFeature feature, int interpolation, bool noinit, bool alltouch) const {
         if (Options::verbose() > 2) std::cout << basename() << " warping into " << imgout.basename() << std::endl;
 
         GeoRaster imgin(*this);
@@ -254,6 +254,8 @@ namespace gip {
             papszOptions = CSLSetNameValue(papszOptions,"INIT_DEST", NULL);
         else
             papszOptions = CSLSetNameValue(papszOptions,"INIT_DEST","NO_DATA");
+        if (alltouch)
+            papszOptions = CSLSetNameValue(papszOptions, "CUTLINE_ALL_TOUCHED", "TRUE");
         papszOptions = CSLSetNameValue(papszOptions,"WRITE_FLUSH","YES");
         papszOptions = CSLSetNameValue(papszOptions,"NUM_THREADS",to_string(Options::cores()).c_str());
         psWarpOptions->papszWarpOptions = papszOptions;
@@ -273,6 +275,7 @@ namespace gip {
             // Create cutline transform to pixel coordinates        
             papszOptionsCutline = CSLSetNameValue( papszOptionsCutline, "DST_SRS", imgout.srs().c_str() );
             papszOptionsCutline = CSLSetNameValue( papszOptionsCutline, "INSERT_CENTER_LONG", "FALSE" );
+
             oTransformer.hSrcImageTransformer = GDALCreateGenImgProjTransformer2( srcDS, NULL, papszOptionsCutline );
             site_t = site->clone();
             site_t->transform(&oTransformer);

GIP/gip/GeoAlgorithms.h

@@ -57,7 +57,8 @@ namespace gip {
     //! Create single image from multiple input images using vector file footprint
     GeoImage cookie_cutter(const std::vector<GeoImage>& geoimgs, std::string filename="",
         GeoFeature feature=GeoFeature(), bool crop=false, std::string proj="",
-        float xres=1.0, float yres=1.0, int interpolation=0, dictionary options=dictionary());
+        float xres=1.0, float yres=1.0, int interpolation=0, dictionary options=dictionary(),
+        bool alltouch=false);
 
     //! Kmeans
     GeoImage kmeans(const GeoImage&, std::string, unsigned int classes=5, unsigned int iterations=5,

GIP/gip/GeoImage.h

@@ -426,7 +426,7 @@ namespace gip {
             bool crop=false, std::string proj="EPSG:4326",
             float xres=1.0, float yres=1.0, int interpolation=0) const;
 
-        GeoImage& warp_into(GeoImage&, GeoFeature=GeoFeature(), int=0, bool=false) const;
+        GeoImage& warp_into(GeoImage&, GeoFeature=GeoFeature(), int=0, bool=false, bool=false) const;
 
     protected:
         //! Vector of raster bands

GIP/gip/GeoRaster.h

@@ -369,7 +369,7 @@ namespace gip {
         template<class T> GeoRaster& write(CImg<T> img, Chunk chunk=Chunk());
         template<class T> GeoRaster& save(GeoRaster& raster) const;
 
-        GeoRaster& warp_into(GeoRaster&, GeoFeature=GeoFeature(), int=0, bool=false) const;
+        GeoRaster& warp_into(GeoRaster&, GeoFeature=GeoFeature(), int=0, bool=false, bool=false) const;
 
         //! NoData mask: 1's where it's bad data
         CImg<unsigned char> nodata_mask(Chunk chunk=Chunk()) const {

setup.py

@@ -38,6 +38,8 @@
 from setuptools.command.build_ext import build_ext
 from setuptools.command.install import install
 from setuptools.command.develop import develop
+from setuptools.command.build_py import build_py
+
 #from wheel.bdist_wheel import bdist_wheel
 from distutils import sysconfig
 
@@ -127,7 +129,7 @@ def run(self):
             os.symlink(os.path.basename(gip_module._file_name), link)
 
         # build the extensionss
-        build_ext.run(self)
+        super().run()
 
         # for mac update runtime library location. Use otool -L to see shared libs in a .so
         if sys.platform == 'darwin':
@@ -149,20 +151,19 @@ def run(self):
 class _develop(develop):
     def run(self):
         log.debug('_develop run')
-        develop.run(self)
+        super().run()
         # move lib files into gippy directory
         [shutil.move(f, 'gippy/') for f in glob.glob('*.so')]
         if sysconfig.get_config_var('SOABI') is not None:
             # rename libgip if Python 3.2+
             os.rename(gip_module._file_name, 'gippy/libgip.so')
 
 
-class _install(install):
+class _build_py(build_py):
     def run(self):
-        log.debug('_install run')
-        # ensure extension built before packaging
-        self.run_command('build_ext')
-        install.run(self)
+        log.debug('_build_py run')
+        self.run_command("build_ext")
+        return super().run()
 
 
 #class _bdist_wheel(bdist_wheel):
@@ -272,7 +273,6 @@ def run(self):
     cmdclass={
         "build_ext": _build_ext,
         "develop": _develop,
-        "install": _install,
-        #"bdist_wheel": _bdist_wheel,
+        "build_py": _build_py,
     }
 )

test/test_GeoAlgorithms.py

@@ -52,36 +52,44 @@ def test_cookiecutter(self):
         res = geoimg1.resolution()
         imgout = alg.cookie_cutter([geoimg1, geoimg2], xres=res.x(), yres=res.y())
         ext = imgout.extent()
-        self.assertEqual(ext.x0(), 0.0)
-        self.assertEqual(ext.y0(), 0.0)
-        self.assertEqual(ext.width(), 2.0)
-        self.assertEqual(ext.height(), 1.0)
+        # This appears to be accurate to 7 decimal places.
+        # Is something getting converted from a double to a float somewhere?
+        self.assertAlmostEqual(ext.x0(), 0.0)
+        self.assertAlmostEqual(ext.y0(), 0.0)
+        self.assertAlmostEqual(ext.width(), 2.0, places=6)
+        self.assertAlmostEqual(ext.height(), 1.0)  # ''
+        self.assertAlmostEqual(imgout.resolution().x(), res.x())
+        self.assertAlmostEqual(imgout.resolution().y(), res.y())
 
     def test_cookiecutter_gain(self):
         """ Cookie cutter on int image with floating point gain """
         bbox = np.array([0.0, 0.0, 1.0, 1.0])
         geoimg = gp.GeoImage.create(xsz=1000, ysz=1000, bbox=bbox, dtype='int16')
         geoimg.set_gain(0.0001)
-        arr = np.zeros((1000,1000)) + 0.0001
+        arr = np.zeros((1000,1000))
         arr[0:500,:] = 0.0002
         geoimg.write(deepcopy(arr))
         res = geoimg.resolution()
         imgout = alg.cookie_cutter([geoimg], xres=res.x(), yres=res.y())
-        np.testing.assert_array_equal(arr, imgout.read())
+        np.testing.assert_array_almost_equal(arr, imgout.read())
 
     def test_cookiecutter_real(self):
         """ Cookie cutter on single real image """
         geoimg = gpt.get_test_image().select(['red']) #, 'green', 'blue'])
+        iext = geoimg.extent()
         vpath = os.path.join(os.path.dirname(__file__), 'vectors')
         # test with feature of different projection
         feature = gp.GeoVector(os.path.join(vpath, 'aoi1_epsg4326.shp'))
         extin = feature.extent()
         imgout = alg.cookie_cutter([geoimg], feature=feature[0], xres=0.0003, yres=0.0003)
         extout = imgout.extent()
-        self.assertAlmostEqual(extout.x0(), extin.x0())
-        self.assertAlmostEqual(extout.y0(), extin.y0())
-        self.assertAlmostEqual(extout.x1(), extin.x1())
-        self.assertAlmostEqual(extout.y1(), extin.y1())
+        self.assertAlmostEqual(extout.x0() + 0.00015, extin.x0())
+        self.assertAlmostEqual(extout.y0() + 0.00015, extin.y0())
+        # cookie cutter will never add more than a pixel and a half in width
+        self.assertTrue(extout.x1() - extin.x1() < 0.0045)
+        self.assertTrue(extout.y1() - extin.y1() < 0.0045)
+        self.assertAlmostEqual(imgout.resolution().x(),  0.0003)
+        self.assertAlmostEqual(imgout.resolution().y(), -0.0003)
 
     def test_cookiecutter_real_reproj(self):
         """ Test with different projection """
@@ -92,10 +100,13 @@ def test_cookiecutter_real_reproj(self):
         # test extent matches feature
         imgout = alg.cookie_cutter([geoimg], feature=feature[0], xres=30.0, yres=30.0)
         extout = imgout.extent()
-        self.assertAlmostEqual(extout.x0(), extin.x0())
-        self.assertAlmostEqual(extout.y0(), extin.y0())
-        self.assertAlmostEqual(extout.x1(), extin.x1())
-        self.assertAlmostEqual(extout.y1(), extin.y1())
+        self.assertAlmostEqual(extout.x0() + 15, extin.x0())
+        self.assertAlmostEqual(extout.y0() + 15, extin.y0())
+        # cookie cutter will never add more than a pixel and a half in width
+        self.assertTrue(extout.x1() - extin.x1() < 45.0)
+        self.assertTrue(extout.y1() - extin.y1() < 45.0)
+        self.assertEqual(imgout.resolution().x(),  30.0)
+        self.assertEqual(imgout.resolution().y(), -30.0)
 
     def test_cookiecutter_real_crop(self):
         """ Test cookie cutter with cropping """
@@ -105,10 +116,11 @@ def test_cookiecutter_real_crop(self):
         imgout = alg.cookie_cutter([geoimg], feature=feature[0], xres=30.0, yres=30.0, crop=True)
         extin = feature.extent()
         extout = imgout.extent()
-        self.assertTrue(extout.x0() >= extin.x0())
-        self.assertTrue(extout.y0() >= extin.y0())
-        self.assertTrue(extout.x1() <= extin.x1())
-        self.assertTrue(extout.y1() <= extin.y1())
+        self.assertTrue(extout.x0() + 15 >= extin.x0())  # half pixel shift
+        self.assertTrue(extout.y0() + 15 >= extin.y0())  # half pixel shift
+        # cookie cutter will never add more than a pixel and a half in width
+        self.assertTrue(extout.x1() - extin.x1() < 45.0)
+        self.assertTrue(extout.y1() - extin.y1() < 45.0)
 
     def test_ndvi(self):
         """ Calculate NDVI using gippy and apply colortable """