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)