Fixes lat/lon extent computation. (9726c9c1) · Commits · aflab / astrogeology / Autocnet

autocnet/fileio/io_gdal.py

+38 −22

Original line number	Diff line number	Diff line
		import json
		import warnings

		import pvl
		@@ -95,6 +94,10 @@ class GeoDataset(object):
		the second value is the lower right corner of the lower right pixel.
		This list is in the form [(minx, miny), (maxx, maxy)].

		xy_corners : list
		of tuple corner coordinates in the form:
		[upper left, lower left, lower right, upper right]

		y_rotation : float
		The geotransform coefficient that represents the rotation about the y-axis.
		Note: This is the fifth value geotransform array.
		@@ -206,15 +209,21 @@ class GeoDataset(object):
		@property
		def latlon_extent(self):
		if not getattr(self, '_latlon_extent', None):
		try:
		if self.footprint:
		fp = self.footprint
		# If we have a footprint, no need to compute pixel to latlon
		# If we have a footprint, do not worry about computing a lat/lon transform
		lowerlat, upperlat, lowerlon, upperlon = fp.GetEnvelope()
		except:
		xy_extent = self.xy_extent
		lowerlat, lowerlon = self.pixel_to_latlon(xy_extent[0][0], xy_extent[0][1])
		upperlat, upperlon = self.pixel_to_latlon(xy_extent[1][0], xy_extent[1][1])
		self._latlon_extent = [(lowerlat, lowerlon), (upperlat, upperlon)]
		self._footprint = [(upperlat, lowerlon),
		(lowerlat, lowerlon),
		(lowerlat, upperlon),
		(upperlat, upperlon)]
		else:
		xy_corners = self.xy_corners
		self._latlon_extent = []
		self.coordinate_transformation
		for x, y in xy_corners:
		x, y = self.pixel_to_latlon(x,y)
		self._latlon_extent.append((x,y))
		return self._latlon_extent

		@property
		@@ -241,20 +250,27 @@ class GeoDataset(object):
		stream = str(f.read(num_polygon_bytes))
		self._footprint = ogr.CreateGeometryFromWkt(stream)
		except:
		# I dislike that this is copied from latlonext, but am unsure
		# how to avoid the cyclical footprint to latlon_extent property hits.
		xy_extent = self.xy_extent
		lowerlat, lowerlon = self.pixel_to_latlon(xy_extent[0][0], xy_extent[0][1])
		upperlat, upperlon = self.pixel_to_latlon(xy_extent[1][0], xy_extent[1][1])
		geom = {"type": "Polygon", "coordinates": [[[lowerlat, lowerlon],
		[lowerlat, upperlon],
		[upperlat, upperlon],
		[upperlat, lowerlon],
		[lowerlat, lowerlon]]]}
		self._footprint = ogr.CreateGeometryFromJson(json.dumps(geom))
		self._footprint = None

		return self._footprint

		@property
		def xy_corners(self):
		if not getattr(self, '_xy_corners', None):
		self._xy_corners = []

		gt = self.geotransform
		x = [0, self.dataset.RasterXSize]
		y = [0, self.dataset.RasterYSize]

		for px in x:
		for py in y:
		xc = gt[0] + (px * gt[1]) + (py * gt[2])
		yc = gt[3] + (px * gt[4]) + (py * gt[5])
		self._xy_corners.append((xc, yc))
		y.reverse()
		return self._xy_corners

		@property
		def xy_extent(self):
		if not getattr(self, '_xy_extent', None):
		@@ -384,9 +400,9 @@ class GeoDataset(object):

		"""
		try:
		geotransform = self.geotransform
		x = geotransform[0] + (x * geotransform[1]) + (y * geotransform[2])
		y = geotransform[3] + (x * geotransform[4]) + (y * geotransform[5])
		#geotransform = self.geotransform
		#x = geotransform[0] + (x * geotransform[1]) + (y * geotransform[2])
		#y = geotransform[3] + (x * geotransform[4]) + (y * geotransform[5])
		lon, lat, _ = self.coordinate_transformation.TransformPoint(x, y)
		except:
		lat = lon = None

autocnet/fileio/tests/test_io_gdal.py

+16 −7

Original line number	Diff line number	Diff line
		@@ -29,7 +29,7 @@ class TestMercator(unittest.TestCase):

		def test_pixel_to_latlon(self):
		lat, lon = self.dataset.pixel_to_latlon(0, 0)
		self.assertAlmostEqual(lat, 55.3322890518, 6)
		self.assertAlmostEqual(lat, 0.0, 6)
		self.assertAlmostEqual(lon, 0.0, 6)

		def test_scale(self):
		@@ -39,8 +39,16 @@ class TestMercator(unittest.TestCase):
		xy_extent = self.dataset.xy_extent
		self.assertEqual(xy_extent, [(0.0, 3921610.0), (10667520.0, -3921610.0)])

		def test_xy_corners(self):
		xy_corners = self.dataset.xy_corners
		self.assertEqual(xy_corners, [(0.0, 3921610.0), (0.0, -3921610.0),
		(10667520.0, -3921610.0), (10667520.0, 3921610.0)])

		def test_latlon_extent(self):
		self.assertEqual(self.dataset.latlon_extent, [(-90, -150.4067721290261), (90.0, 0.0)])
		self.assertEqual(self.dataset.latlon_extent, [(55.33228905180849, 0.0),
		(-55.3322890518085, 0.0),
		(-55.3322890518085, 179.96751473604124),
		(55.33228905180849, 179.96751473604124)])

		def test_spheroid(self):
		sphere = self.dataset.spheroid
		@@ -85,6 +93,7 @@ class TestMercator(unittest.TestCase):
		self.assertEqual(arr.dtype, np.int8)
		self.assertAlmostEqual(np.mean(arr), 10.10353227, 6)


		class TestLambert(unittest.TestCase):
		def setUp(self):
		self.dataset = io_gdal.GeoDataset(get_path('Lunar_LRO_LOLA_Shade_MAP2_90.0N20.0_LAMB.tif'))
		@@ -104,8 +113,8 @@ class TestLambert(unittest.TestCase):

		def test_pixel_to_latlon(self):
		lat, lon = self.dataset.pixel_to_latlon(0,0)
		self.assertAlmostEqual(lat, 69.90349154912009, 6)
		self.assertAlmostEqual(lon, -29.72166902463681, 6)
		self.assertAlmostEqual(lat, 90.0, 6)
		self.assertAlmostEqual(lon, 20.0, 6)

		def test_latlon_to_pixel(self):
		lat, lon = 69.90349154912009, -29.72166902463681
		@@ -141,8 +150,8 @@ class TestPolar(unittest.TestCase):

		def test_pixel_to_latlon(self):
		lat, lon = self.dataset.pixel_to_latlon(0,0)
		self.assertAlmostEqual(lat, 42.2574735013, 6)
		self.assertAlmostEqual(lon, -135.0, 6)
		self.assertAlmostEqual(lat, 90, 6)
		self.assertAlmostEqual(lon, 0.0, 6)

		def test_latlon_to_pixel(self):
		lat, lon = 42.2574735013, -135.0