diff --git a/.circleci/config.yml b/.circleci/config.yml index 45536cd..7ce7423 100644 --- a/.circleci/config.yml +++ b/.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 diff --git a/.gitignore b/.gitignore index 1532d21..3934d4c 100644 --- a/.gitignore +++ b/.gitignore @@ -14,3 +14,4 @@ gippy/gippy_wrap.cpp gippy/algorithms.py gippy/algorithms_wrap.cpp /.vs +.venv* diff --git a/Dockerfile b/Dockerfile index 010304f..bf4035b 100644 --- a/Dockerfile +++ b/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 diff --git a/GIP/GeoAlgorithms.cpp b/GIP/GeoAlgorithms.cpp index 3d9be1a..eef306b 100644 --- a/GIP/GeoAlgorithms.cpp +++ b/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& 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,21 @@ 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 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; + + // shift lower left corner IFF transition from vector to raster space + 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 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 +288,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 2) std::cout << basename() << " warping into " << imgout.basename() << std::endl; // Assume that these have the same number of bands for (unsigned int b=0; bnbands(); b++) { - (*this)[b].warp_into(imgout[b], feature, interpolation, noinit); + (*this)[b].warp_into(imgout[b], feature, interpolation, noinit, alltouch); } return imgout; diff --git a/GIP/GeoRaster.cpp b/GIP/GeoRaster.cpp index 75e9356..5331334 100644 --- a/GIP/GeoRaster.cpp +++ b/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); diff --git a/GIP/gip/GeoAlgorithms.h b/GIP/gip/GeoAlgorithms.h index 5b7cf56..a297c82 100644 --- a/GIP/gip/GeoAlgorithms.h +++ b/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& 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, diff --git a/GIP/gip/GeoImage.h b/GIP/gip/GeoImage.h index 88c445a..07765b2 100644 --- a/GIP/gip/GeoImage.h +++ b/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 diff --git a/GIP/gip/GeoRaster.h b/GIP/gip/GeoRaster.h index c31895b..67ba604 100644 --- a/GIP/gip/GeoRaster.h +++ b/GIP/gip/GeoRaster.h @@ -369,7 +369,7 @@ namespace gip { template GeoRaster& write(CImg img, Chunk chunk=Chunk()); template 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 nodata_mask(Chunk chunk=Chunk()) const { diff --git a/setup.py b/setup.py index fed80c0..a5a73f2 100755 --- a/setup.py +++ b/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,7 +151,7 @@ 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: @@ -157,12 +159,11 @@ def run(self): 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, } ) diff --git a/test/test_GeoAlgorithms.py b/test/test_GeoAlgorithms.py index 52b78f4..34becd9 100644 --- a/test/test_GeoAlgorithms.py +++ b/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 """