Biomass reproject + update links#575
Conversation
|
Check out this pull request on See visual diffs & provide feedback on Jupyter Notebooks. Powered by ReviewNB |
hrodmn
left a comment
There was a problem hiding this comment.
@HarshiniGirish I'm sorry this PR slipped past me a few months ago!
The STAC item bbox approach is clever but we don't need to do that! I am not sure why these files have a null crs but they are georeferenced. Instead of a crs they have gcps defined (gcps = Ground Control Points).
Here is how you can use it to do the reprojection without guessing using the STAC metadata:
rio_env = rasterio.Env(GDAL_HTTP_HEADERS=f"Authorization: Bearer {token}")
with rio_env:
with rasterio.open(asset_href) as ds:
arr = ds.read(1)
gcps_list, gcp_crs = ds.gcps
valid = arr[np.isfinite(arr)]
print("Shape:", arr.shape)
print("CRS:", ds.crs)
print("GCPS list:", gcps_list)
print("GCP CRS:", gcp_crs)
print("Transform:", ds.transform)
print("Min/Max:", float(valid.min()), float(valid.max()))
src_arr = arr.copy()
src_height, src_width = src_arr.shape
# Reproject to EPSG:4326 on a clean output grid
dst_crs = "EPSG:4326" # don't use EPSG:4326 for raster data if possible!
dst_transform, dst_width, dst_height = calculate_default_transform(
src_crs=gcp_crs,
dst_crs=dst_crs,
width=src_width,
height=src_height,
gcps=gcps_list
)
reprojected = np.empty((dst_height, dst_width), dtype=src_arr.dtype)
reproject(
source=src_arr,
destination=reprojected,
gcps=gcps_list,
src_crs=gcp_crs,
dst_transform=dst_transform,
dst_crs=dst_crs,
resampling=Resampling.nearest
)
print("Destination CRS:", dst_crs)
print("Reprojected shape:", reprojected.shape)
print("Destination transform:", dst_transform)I think we should encourage users to use a different CRS than epsg:4326 for array data like this, so you could figure out roughly where the sample file is then pick a projected CRS that fits the area.
|
thanks for the feedback @hrodmn. It's up for review again! |
No description provided.