Commit fdc06e7c authored by acpaquette's avatar acpaquette Committed by jlaura
Browse files

North_up/footprint (#50)

* Removed north_up and added basic footprint creation.

* Adjusted geotransform calculation to get a geo before trying to compute it.

* Added tests for latlon_corners and extent calculation.

* potential travis fix

* Other potential fix.
parent cf4eb7d2
Loading
Loading
Loading
Loading
+1 −0
Original line number Diff line number Diff line
@@ -4,6 +4,7 @@ sudo: false
branches:
  only:
    - master

os:
  - linux
  - osx
+30 −21
Original line number Diff line number Diff line
@@ -184,7 +184,9 @@ class GeoDataset(object):
        top left y, y-rotation, n-s pixel resolution]
        """
        if not getattr(self, '_geotransform', None):
            if self.footprint:
            try:
                self._geotransform = self.dataset.GetGeoTransform()
            except:
                coords = json.loads(self.footprint.ExportToJson())['coordinates'][0][0]
                ul, ur, lr, ll = geofuncs.find_four_corners(coords)

@@ -196,8 +198,7 @@ class GeoDataset(object):
                yskew = (ur[1] - ystart) / xsize
                yscale = (ll[1] - ystart) / ysize
                self._geotransform = [xstart, xscale, xskew, ystart, yskew, yscale]
            else:
                self._geotransform = self.dataset.GetGeoTransform()

        return self._geotransform

    @property
@@ -225,13 +226,6 @@ class GeoDataset(object):
            self._unit_type = self.dataset.GetRasterBand(1).GetUnitType()
        return self._unit_type

    @property
    def north_up(self):
        if self.footprint:
            return geofuncs.is_clockwise(json.loads(self.footprint.ExportToJson())['coordinates'][0][0])
        else:
            return True

    @property
    def spatial_reference(self):
        if not getattr(self, '_srs', None):
@@ -289,6 +283,25 @@ class GeoDataset(object):
            except:
                self._footprint = None

            try:
                # Get the lat lon corners
                lat = [i[0] for i in self.latlon_corners]
                lon = [i[1] for i in self.latlon_corners]

                # Compute a ogr geometry for the tiff which
                # provides leverage for overlaps
                ring = ogr.Geometry(ogr.wkbLinearRing)
                for point in [*zip(lon, lat)]:
                    ring.AddPoint(*point)
                ring.AddPoint(lon[0], lat[0])

                poly = ogr.Geometry(ogr.wkbPolygon)
                poly.AddGeometry(ring)
                poly.FlattenTo2D()
                self._footprint = poly
            except:
                self._footprint = None

        return self._footprint

    @property
@@ -313,11 +326,13 @@ class GeoDataset(object):
                self._latlon_extent = [(minx, miny),
                                       (maxx, maxy)]
            else:
                # Attempt to compute a basic bounding box footprint
                self._latlon_extent = []
                for x, y in self.xy_extent:
                    x, y = self.pixel_to_latlon(x,y)

                    self._latlon_extent.append((x,y))
                xy_extent = self.xy_extent
                maxlon, minlat = self.pixel_to_latlon(*xy_extent[0])
                minlon, maxlat = self.pixel_to_latlon(*xy_extent[1])
                self._latlon_extent.append((minlat, minlon))
                self._latlon_extent.append((maxlat, maxlon))
        return self._latlon_extent

    @property
@@ -496,15 +511,10 @@ class GeoDataset(object):

        if not pixels:
            array = band.ReadAsArray().astype(dtype)
            if self.north_up == False:
                array = np.flipud(array)
        else:
            # Check that the read start is not outside of the image
            xstart, ystart, xcount, ycount = pixels
            xmax, ymax = map(int, self.xy_extent[1])
            # If the image is south up, flip the roi
            if self.north_up == False:
                ystart = ymax - (ystart + ycount)
            if xstart < 0:
                xstart = 0

@@ -517,8 +527,7 @@ class GeoDataset(object):
            if ystart + ycount > ymax:
                ycount = ymax - ystart
            array = band.ReadAsArray(xstart, ystart, xcount, ycount).astype(dtype)
            #if self.north_up == False:
            #    array = np.flipud(array)

        return array

    def compute_overlap(self, geodata, **kwargs):
+14 −1
Original line number Diff line number Diff line
@@ -50,7 +50,12 @@ class TestMercator(unittest.TestCase):
    def test_latlon_extent(self):
        self.assertEqual(self.dataset.latlon_extent, [(55.33228905180849, 0.0),
                                                      (-55.3322890518085, 179.96751473604124)])

    def test_latlon_corners(self):
        self.assertEqual(self.dataset.latlon_corners, [(55.33228905180849, 0.0),
                                                      (-55.3322890518085, 179.96751473604124)])
    """

    def test_spheroid(self):
        sphere = self.dataset.spheroid
        self.assertAlmostEqual(sphere[0], 3396190.0, 6)
@@ -111,6 +116,14 @@ class TestLambert(unittest.TestCase):
    def test_get_xy_extent(self):
        self.assertEqual(self.dataset.xy_extent, [(0, 0), (239, 275)])

    def test_latlon_extent(self):
        self.assertEqual(self.dataset.latlon_extent, [(-29.721669024636785, 37.834604457982444),
                                                      (69.44231975603968, 69.99086737565507)])
    def test_latlon_corners(self):
        self.assertEqual(self.dataset.latlon_corners, [(-29.721669024636785, 69.99086737565507),
                                                       (-29.721669024636785, 37.834604457982444),
                                                       (69.44231975603968, 37.834604457982444),
                                                       (69.44231975603968, 69.99086737565507)])
    def test_get_no_data_value(self):
        self.assertEqual(self.dataset.no_data_value, 0.0)