Tweak antimeridian hack and timezone handling.

datacube_ows/index/postgis/api.py

@@ -65,10 +65,13 @@ def _query(self,
                ) -> dict[str, Any]:
         query: dict[str, Any] = {}
         if geom:
-            geopoly = self._prep_geom(layer, geom).to_crs("epsg:4326")
-            geopoly = Geometry(fix_shape(geopoly.geom), crs="epsg:4326")
-            query["geopolygon"] = geopoly
-            _LOG.warning("geopolygon is %s", repr(query["geopolygon"]))
+            if geom.crs in layer.dc.index.spatial_indexes():
+                query["geopolygon"] = geom
+            else:
+                # Default to 4326 and take a long hard look at yourself.
+                geopoly = self._prep_geom(layer, geom).to_crs("epsg:4326")
+                geopoly = Geometry(fix_shape(geopoly.geom), crs="epsg:4326")
+                query["geopolygon"] = geopoly
         if products is not None:
             query["product"] = [p.name for p in products]
         if times is not None:

datacube_ows/time_utils.py

@@ -112,7 +112,7 @@ def local_solar_date_range(geobox: GeoBox, date: datetime.date) -> tuple[datetim
     :param date: A date object
     :return: A tuple of two UTC datetime objects, spanning 1 second shy of 24 hours.
     """
-    tz: datetime.tzinfo = tz_for_geometry(geobox.geographic_extent)
+    tz: datetime.tzinfo = tz_for_geometry(geobox.extent)
     start = datetime.datetime(date.year, date.month, date.day, 0, 0, 0, tzinfo=tz)
     end = datetime.datetime(date.year, date.month, date.day, 23, 59, 59, tzinfo=tz)
     return (start.astimezone(utc), end.astimezone(utc))
@@ -174,20 +174,14 @@ def tz_for_geometry(geom: Geometry) -> datetime.tzinfo:
     :return: A timezone object
     """
     crs_geo = CRS("EPSG:4326")
-    geo_geom: Geometry = geom.to_crs(crs_geo)
-    centroid: Geometry = geo_geom.centroid
+    raw_centroid = geom.centroid
+    centroid: Geometry = raw_centroid.to_crs(crs_geo)
     try:
         # 1. Try being smart with the centroid of the geometry
         return tz_for_coord(centroid.coords[0][0], centroid.coords[0][1])
     except NoTimezoneException:
         pass
-    for pt in geo_geom.boundary.coords:
-        try:
-            # 2. Try being smart all the points in the geometry
-            return tz_for_coord(pt[0], pt[1])
-        except NoTimezoneException:
-            pass
-    # 3. Meh, just use longitude
+    # 2. Meh, just use longitude
     offset = round(centroid.coords[0][0] / 15.0)
     return datetime.timezone(datetime.timedelta(hours=offset))
 

datacube_ows/wms_utils.py

@@ -78,21 +78,14 @@ def _get_geobox(args, crs):
     if minx == maxx or miny == maxy:
         raise WMSException("Bounding box must enclose a non-zero area")
 
-    if crs.epsg == 3857:
-        # Web Mercator anti-meridian hack
-        if maxx < -13_000_000 or minx > 13_000_000:
-            _LOG.warning("Applying anti-meridian hack! x=~ %f, %f    y=~ %f, %f", minx, maxx, miny, maxy)
-            # Closer to the anti-meridian than the prime meridian:
+    if crs.epsg == 3857 and (maxx < -13_000_000 or minx > 13_000_000):
+            # EPSG:3857 query AND closer to the anti-meridian than the prime meridian:
             # re-project to epsg:3832 (Pacific Web-Mercator)
             ll = geom.point(x=minx, y=miny, crs=crs).to_crs("epsg:3832")
             ur = geom.point(x=maxx, y=maxy, crs=crs).to_crs("epsg:3832")
             minx, miny = ll.coords[0]
             maxx, maxy = ur.coords[0]
             crs = geom.CRS("epsg:3832")
-        else:
-            _LOG.warning("NOT applying anti-meridian hack! x=~ %f, %f  y=~ %f %f", minx, maxx, miny, maxy)
-    else:
-        _LOG.warning("Not a 3857 query!")
 
     return create_geobox(
         crs,
@@ -118,11 +111,11 @@ def zoom_factor(args, crs):
     # Project to a geographic coordinate system
     # This is why we can't just use the regular geobox.  The scale needs to be
     # "standardised" in some sense, not dependent on the CRS of the request.
-    geo_crs = geom.CRS("EPSG:4326")
+    # TODO: can we do better in polar regions?
     minx, miny, maxx, maxy = _bounding_pts(
         minx, miny,
         maxx, maxy,
-        crs, dst_crs=geo_crs
+        crs, dst_crs="epsg:4326"
     )
     # Create geobox affine transformation (N.B. Don't need an actual Geobox)
     affine = Affine.translation(minx, miny) * Affine.scale((maxx - minx) / width, (maxy - miny) / height)
@@ -508,8 +501,7 @@ def solar_correct_data(data, dataset):
     native_x = (dataset.bounds.right + dataset.bounds.left) / 2.0
     native_y = (dataset.bounds.top + dataset.bounds.bottom) / 2.0
     pt = geom.point(native_x, native_y, dataset.crs)
-    crs_geo = geom.CRS("EPSG:4326")
-    geo_pt = pt.to_crs(crs_geo)
+    geo_pt = pt.to_crs("epsg:4326")
     data_time = dataset.center_time.astimezone(utc)
     data_lon, data_lat = geo_pt.coords[0]
 

setup.py

@@ -10,7 +10,7 @@
 from setuptools import find_packages, setup
 
 install_requirements = [
-    'datacube[performance,s3]>=1.9.0-rc9',
+    'datacube[performance,s3]>=1.9.0-rc11',
     'flask',
     'requests',
     'affine',