Merge pull request #2 from TUW-GEO/develop

merge develop with master

.travis.yml

@@ -0,0 +1,48 @@
+language: python
+sudo: false
+addons:
+  apt:
+    packages:
+      - gfortran
+      - gcc
+      - libgrib-api-dev
+notifications:
+  email: false
+python:
+  # We don't actually use the Travis Python, but this keeps it organized.
+  - "2.7"
+  - "3.4"
+  - "3.5"
+  - "3.6"
+install:
+  # You may want to periodically update this, although the conda update
+  # conda line below will keep everything up-to-date.  We do this
+  # conditionally because it saves us some downloading if the version is
+  # the same.
+  - if [[ "$TRAVIS_PYTHON_VERSION" == "2.7" ]]; then
+      wget http://repo.continuum.io/miniconda/Miniconda-latest-Linux-x86_64.sh -O miniconda.sh;
+    else
+      wget http://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh;
+    fi
+  - bash miniconda.sh -b -p $HOME/miniconda
+  - export PATH="$HOME/miniconda/bin:$PATH"
+  - hash -r
+  - conda config --set always_yes yes --set changeps1 no
+  - conda update -q conda
+  # Useful for debugging any issues with conda
+  - conda info -a
+
+  - conda create -q -n test-environment -c conda-forge python=$TRAVIS_PYTHON_VERSION numpy scipy pytest pip gdal
+  pyproj numba astropy toolz dask
+  - source activate test-environment
+  - python setup.py install
+  - pip list
+  - which pip
+  - which python
+
+script:
+  - python setup.py test
+after_success:
+  # report coverage results to coveralls.io
+  - pip install coveralls
+  - coveralls

conda_environment.yml

@@ -0,0 +1,49 @@
+name: C:\PythonMiniconda27_4.3.21\envs\TPS
+channels: !!python/tuple
+- !!python/unicode
+  'defaults'
+dependencies:
+- astropy=2.0.1=np113py27_0
+- certifi=2016.2.28=py27_0
+- colorama=0.3.9=py27_0
+- coverage=4.4.1=py27_0
+- curl=7.52.1=vc9_0
+- enum34=1.1.6=py27_0
+- funcsigs=1.0.2=py27_0
+- gdal=2.1.0=py27_0
+- geos=3.5.0=vc9_0
+- hdf4=4.2.12=vc9_1
+- hdf5=1.8.15.1=vc9_4
+- jpeg=9b=vc9_0
+- kealib=1.4.6=vc9_0
+- libgdal=2.1.0=vc9_0
+- libnetcdf=4.3.3.1=vc9_4
+- libtiff=4.0.6=vc9_3
+- llvmlite=0.19.0=py27_0
+- mkl=2017.0.3=0
+- numba=0.34.0=np113py27_0
+- numpy=1.13.1=py27_0
+- pip=9.0.1=py27_1
+- proj4=4.9.2=vc9_0
+- py=1.4.34=py27_0
+- pytest=3.2.1=py27_0
+- pytest-cov=2.3.1=py27_0
+- python=2.7.13=1
+- scipy=0.19.1=np113py27_0
+- setuptools=36.4.0=py27_0
+- singledispatch=3.4.0.3=py27_0
+- six=1.10.0=py27_0
+- vs2008_runtime=9.00.30729.5054=0
+- vs2015_runtime=14.0.25420=0
+- wheel=0.29.0=py27_0
+- wincertstore=0.2=py27_0
+- xerces-c=3.1.4=vc9_0
+- zlib=1.2.11=vc9_0
+- pip:
+  - click==6.7
+  - dask==0.15.3
+  - equi7grid (c:\code\tps\equi7grid)==0.0.post0.dev1+ng28b01b3.dirty
+  - pyproj==1.9.5.1
+  - pytileproj (c:\code\tps\pytileproj)==0.0.post0.dev10+nfe209d0.dirty
+  - toolz==0.8.2
+prefix: C:\PythonMiniconda27_4.3.21\envs\TPS

pytileproj/TiledProjection.py → pytileproj/base.py

@@ -38,10 +38,10 @@
 '''
 
 import abc
+import itertools
 
 import numpy as np
 import dask.array as da
-
 from osgeo import osr
 
 import geometry as geometry
@@ -203,7 +203,7 @@ def locate_geometry_in_subgrids(self, geom):
         return covering_subgrid
 
 
-    def lonlat2xy(self, lat, lon, subgrid=None):
+    def lonlat2xy(self, lon, lat, subgrid=None):
         '''
         converts latitude and longitude coordinates to TPS grid coordinates
 
@@ -216,13 +216,13 @@ def lonlat2xy(self, lat, lon, subgrid=None):
         # TODO use pyproj.transform for performance when having many points
         if subgrid is None:
             vfunc = np.vectorize(self._lonlat2xy)
-            return vfunc(lat, lon)
+            return vfunc(lon, lat)
         else:
             vfunc = np.vectorize(self._lonlat2xy_subgrid)
-            return vfunc(lat, lon, subgrid)
+            return vfunc(lon, lat, subgrid)
 
 
-    def _lonlat2xy(self, lat, lon):
+    def _lonlat2xy(self, lon, lat):
         """
         finds overlapping subgrids of given geometry and computes the projected coordinates
         """
@@ -241,7 +241,7 @@ def _lonlat2xy(self, lat, lon):
         return np.full_like(x, subgrid, dtype=(np.str, len(subgrid))), x, y
 
 
-    def _lonlat2xy_subgrid(self, lat, lon, subgrid):
+    def _lonlat2xy_subgrid(self, lon, lat, subgrid):
         '''
         computes the projected coordinates in given subgrid
 
@@ -272,6 +272,143 @@ def get_tiletype(self, res):
     def get_tilesize(self, res):
         pass
 
+    def search_tiles_in_geo_roi(self,
+                                geom_area=None,
+                                extent=None,
+                                epsg=4326,
+                                wkt=None,
+                                subgrid_ids=None,
+                                coverland=False,
+                                gdal_path=None):
+
+        """
+        Search the tiles which are intersected by the poly_roi area.
+
+        Parameters
+        ----------
+        geom_area : geometry
+            a polygon or multipolygon geometery object representing the ROI
+        extent : list
+            It is a polygon representing the rectangle region of interest
+            in the format of [xmin, ymin, xmax, ymax].
+        epsg : str
+            EPSG CODE defining the spatial reference system, in which
+            the geometry or extent is given. Default is LatLon (EPSG:4326)
+        subgrid_ids : list
+            grid ID to specified which continents you want to search. Default
+            value is None for searching all continents.
+
+        Returns
+        -------
+        list
+            return a list of  the overlapped tiles' name.
+            If not found, return empty list.
+        """
+        # check input grids
+        if subgrid_ids is None:
+            subgrid_ids = self.subgrids.keys()
+        else:
+            subgrid_ids = [x.upper() for x in subgrid_ids]
+            if set(subgrid_ids).issubset(set(self.subgrids.keys())):
+                subgrid_ids = list(subgrid_ids)
+            else:
+                raise ValueError("Invalid agrument: grid must one of [ %s ]." %
+                                 " ".join(self.subgrids.keys()))
+
+        if not geom_area and not extent:
+            print "Error: either geom or extent should be given as the ROI."
+            return list()
+
+        # obtain the geometry of ROI
+        if not geom_area:
+            if epsg is not None:
+                geom_area = geometry.extent2polygon(extent, epsg=epsg)
+            elif wkt is not None:
+                geom_area = geometry.extent2polygon(extent, wkt=wkt)
+            else:
+                raise ValueError("Error: either epsg or wkt of ROI coordinates must be defined!")
+
+        # load lat-lon spatial reference as the default
+        geo_sr = TPSProjection(epsg=4326).osr_spref
+
+        geom_sr = geom_area.GetSpatialReference()
+        if geom_sr is None:
+            geom_area.AssignSpatialReference(geo_sr)
+        elif not geom_sr.IsSame(geo_sr):
+            geom_area.TransformTo(geo_sr)
+
+        # intersect the given grid ids and the overlapped ids
+        overlapped_grids = self.locate_geometry_in_subgrids(geom_area)
+        subgrid_ids = list(set(subgrid_ids) & set(overlapped_grids))
+
+        # finding tiles
+        overlapped_tiles = list()
+        for sgrid_id in subgrid_ids:
+            overlapped_tiles.extend(
+                self._search_sgrid_tiles(geom_area, sgrid_id, coverland))
+        return overlapped_tiles
+
+    def _search_sgrid_tiles(self, geom, sgrid_id, coverland):
+        """
+        Search the tiles which are overlapping with the given grid.
+
+        Parameters
+        ----------
+        area_geometry : geometry
+            It is a polygon geometry representing the region of interest.
+        sgrid_id : string
+            sub grid ID to specified which continent you want to search.
+            Default value is None for searching all continents.
+
+        Returns
+        -------
+        list
+            Return a list of  the overlapped tiles' name.
+            If not found, return empty list.
+        """
+
+        subgrid = getattr(self, sgrid_id)
+
+        # get the spatial reference of the subgrid
+        grid_sr = subgrid.projection.osr_spref
+
+        # get the intersection of the area of interest and grid zone
+        geom.TransformTo(grid_sr)
+
+        intersect = geom.Intersection(geom.Intersection(subgrid.polygon_proj))
+        if not intersect:
+            return list()
+        # The spatial reference need to be set again after intersection
+        #intersect.AssignSpatialReference(geom.GetSpatialReference())
+
+        # transform the area of interest to the grid coordinate system
+        #grid_sr = getattr(self, sgrid_id).projection.osr_spref
+        #intersect.TransformTo(grid_sr)
+
+        # get envelope of the Geometry and cal the bounding tile of the
+        envelope = intersect.GetEnvelope()
+        x_min = int(envelope[0]) / self.core.tile_xsize_m * self.core.tile_xsize_m
+        x_max = (int(envelope[1]) / self.core.tile_xsize_m + 1) * self.core.tile_xsize_m
+        y_min = int(envelope[2]) / self.core.tile_ysize_m * self.core.tile_ysize_m
+        y_max = (int(envelope[3]) / self.core.tile_ysize_m + 1) * self.core.tile_ysize_m
+
+        # make sure x_min and y_min greater or equal 0
+        x_min = 0 if x_min < 0 else x_min
+        y_min = 0 if y_min < 0 else y_min
+
+        # get overlapped tiles
+        overlapped_tiles = list()
+        for x, y in itertools.product(xrange(x_min, x_max, self.core.tile_xsize_m),
+                                      xrange(y_min, y_max, self.core.tile_ysize_m)):
+            geom_tile = geometry.extent2polygon((x, y, x + self.core.tile_xsize_m,
+                                                 y + self.core.tile_xsize_m))
+            if geom_tile.Intersects(intersect):
+                ftile = subgrid.tilesys.point2tilename(x, y)
+                if not coverland or subgrid.tilesys.check_land_coverage(ftile):
+                    overlapped_tiles.append(ftile)
+
+        return overlapped_tiles
+
 class TiledProjection(object):
     """
     Class holding the projection and tiling definition of a tiled projection space.
@@ -406,7 +543,7 @@ def decode_tilename(self, name):
         return a
 
     @abc.abstractmethod
-    def identify_tiles_overlapping_bbox(self, bbox):
+    def identify_tiles_overlapping_xybbox(self, bbox):
         """Light-weight routine that returns
            the name of tiles intersecting the bounding box.
 
@@ -434,7 +571,7 @@ def identify_tiles_overlapping_bbox(self, bbox):
         return tilenames
 
 
-    def create_tiles_overlapping_bbox(self, bbox):
+    def create_tiles_overlapping_xybbox(self, bbox):
         """Light-weight routine that returns
            the name of tiles intersecting the bounding box.
 
@@ -451,7 +588,7 @@ def create_tiles_overlapping_bbox(self, bbox):
             with .subset() not exceeding the bounding box.
 
         """
-        tilenames = self.identify_tiles_overlapping_bbox(bbox)
+        tilenames = self.identify_tiles_overlapping_xybbox(bbox)
         tiles = list()
 
         for t in tilenames:
@@ -479,9 +616,9 @@ def create_tiles_overlapping_bbox(self, bbox):
         return tiles
 
 
-    def create_daskarray_overlapping_bbox(self, bbox):
+    def create_daskarray_overlapping_xybbox(self, bbox):
 
-        tiles = self.create_tiles_per_bbox(bbox)
+        tiles = self.create_tiles_overlapping_xybbox(bbox)
         tilenames = [x.name for x in tiles]
 
         x_anchors = set([t.llx for t in tiles])

pytileproj/PixelDownsampler.py → pytileproj/downsample.py

@@ -49,7 +49,8 @@
 from scipy import ndimage
 from numba import jit
 from astropy.convolution import convolve
-from Kernel import Kernel
+
+from pytileproj.kernel import Kernel
 
 
 class PixelDownsampler(object):
@@ -262,11 +263,11 @@ def calc_pixel_index_pattern(spacing_fine, spacing_coarse):
     return pattern
 
 
-def translate_pixelmaps(spacing_fine, spacing_coarse, target_bbox):
+def translate_pixelmaps(spacing_fine, spacing_coarse, target_bbox, correct_boundary=False):
 
     #TODO:
     # implement ghosting-reader (read buffer around target area), then set this to True
-    correct_boundary = False
+    # correct_boundary = True
 
     # resolution of grid
     res_f = spacing_fine