0

I face a crucial problem with ODC. I try to do a spatial join with a shapefile (polygon drill). I have followed the following steps:

  • Read the shapefile and get the first polygon to test
  • Ensure that projection of the polygon is the same with the polygon of the raw data indexed in the ODC
  • Make a query based on the geopolygon attribute.

The query does return an empty xarray. However, when I remove the geopolygon attribute from the query it does return all the data. Any idea?

import datacube
import math
# import folium
import calendar
import ipywidgets
import numpy as np
# import geopandas as gpd
import matplotlib as mpl
import matplotlib.cm as cm
from matplotlib import colors as mcolours
import matplotlib.patheffects as PathEffects
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from datetime import datetime
from pyproj import Proj, transform
from datacube.utils import geometry
from datacube.utils.geometry import CRS,Geometry
import shapely
from shapely.geometry import shape
import fiona
import rasterio.features
from fiona.crs import from_epsg
dc = datacube.Datacube(app="test",config="/home/noa/datacube.conf")
shape_file = '/home/noa/Desktop/cyprus36SWDparcels3857.shp'
ds = fiona.open(shape_file)
crs = geometry.CRS(ds.crs_wkt)
for f in ds:
    feature_geom = f['geometry']
    print(feature_geom['extent'])
    geom = geometry.Geometry(geom=feature_geom,crs=crs)
    bounds = shape(feature_geom).bounds
    if 'MULTIPOLYGON' in geom.wkt:
        continue
    break
product = 'sentinel1' # 's2a_sen2cor_granule'
query = {
        'geopolygon': geom,
        'time': ('2018-01-01', '2018-01-31')
}
data = dc.load(product=product, output_crs="EPSG:3857",measurements=['vv'],resolution=(-10,10),**query)
Improve this question
3
  • Are you sure that the first geometry is the one you really want? It's not uncommon to see multiple geometries in a shapefile. You can do some debugging by either dumping the contents of the first_geometry or by building one that you know is correct (as in the answer below)
    – A C
    Mar 17, 2021 at 0:13
  • Thanks for the reply. This is just an exampe. I do want all the geometries but I tried only for the first one and a series of random geometries of the shapefile. Then I tried to test all of them and it does not work with any of them.. Note here that the shapefile is within the bounds of the sentinel images indexed in the cube.
    – Thanassis Drivas
    Mar 21, 2021 at 12:27
  • In addition, debug took place: polygons' extent is corrent, images' extent is correct. I tried also with a Polygon that I created from scratch by defining certain coordinates and it does not work. Moreover, I tried changing multiple crs to the Polygons' coordinates. Not working. I am desperate... General Question: Is the coordinates of the asked polygon transformed during the query so to match the projection of the indexed image or must we set the project's coordinate system to the one of the image ?? (FYI I tried also the second one, no luck)
    – Thanassis Drivas
    Mar 21, 2021 at 12:35

1 Answer 1

Reset to default
0

I'd check your CRS. Here's a working query using geopolygon:

import datacube
from datacube.utils.geometry import Geometry, CRS

# Create a Datacube object
dc = datacube.Datacube()

# Wrangle some geometry
geojson = {'type': 'Polygon',
  'coordinates': [[[2304000.0, -3840000.0],
    [2304000.0, -3936000.0],
    [2400000.0, -3936000.0],
    [2400000.0, -3840000.0],
    [2304000.0, -3840000.0]]]}
crs = CRS("epsg:6933")
geom = Geometry(geom=geojson, crs=crs)

# Do a DC Load
dc.load(product=product, time="2020-01-01", geopolygon=geom, output_crs="epsg:4326", resolution=(0.001, 0.001))

This outputs an expected dataset:

Dimensions: latitude: 886 longitude: 995 time: 1

Improve this answer
1
  • Thanks for the reply. However, I have checked both CRS and geometry. There are fine. In addition, there is an unusual behaviour; when I set the a pair of coordinates in the X, Y parameters in the query it does work. But when I remove this coords and pass the geopolygon, it does give an empty array. I have to add here that the product in the ODC database does not carry projection. I dont think that is a problem as I make spatial queries with the aforementioned way. To be safe, I updated the json field in the db by adding the projection, but the "geopolygon" problem remains parameter.
    – Thanassis Drivas
    Mar 11, 2021 at 10:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.