EOA-team / eodal

Earth Observation Data Analysis Library
https://eodal.readthedocs.io/en/latest/
GNU General Public License v3.0
93 stars 15 forks source link

numpy masked array arithmetics; ValueError: output array is read-only for numpy>=1.22 #3

Closed m-i-g-u-e-l closed 2 years ago

m-i-g-u-e-l commented 2 years ago

File "/home/miguel/code/agroscope/eodal/eodal-venv39_v2/lib/python3.9/site-packages/numpy/lib/function_base.py", line 4489, in _lerp subtract(b, diff_b_a * (1 - t), out=lerp_interpolation, where=t >= 0.5)

ValueError: output array is read-only

lerp_interpolation Out [3]: masked

type(lerp_interpolation) Out [4]: <class 'numpy.ma.core.MaskedConstant'>

Example code:

!/usr/bin/env python3

-- coding: utf-8 --

""" Created on Tue Aug 9 11:53:46 2022

@author: miguel """

import cv2

import numpy as np import os

from pathlib import Path from eodal.core.spectral_indices import SpectralIndices from eodal.core.sensors import Sentinel2

make plots larger by default

import matplotlib.pyplot as plt plt.rcParams['figure.figsize'] = [10, 10]

download test data (if not done yet)

import requests from eodal.downloader.sentinel2.utils import unzip_datasets

URL to the public dataset

url = 'https://data.mendeley.com/public-files/datasets/ckcxh6jskz/files/e97b9543-b8d8-436e-b967-7e64fe7be62c/file_downloaded'

testdata_dir = Path('../../data') if not os.path.exists(testdata_dir): os.makedirs(testdata_dir) testdata_fname = testdata_dir.joinpath('S2A_MSIL2A_20190524T101031_N0212_R022_T32UPU_20190524T130304.zip') testdata_fname_unzipped = Path(testdata_fname.as_posix().replace('.zip', '.SAFE'))

check first if the dataset has been already downloaded; only start the download if the dataset is not yet available locally

if not testdata_fname_unzipped.exists():

# download dataset
r = requests.get(url, stream=True)
r.raise_for_status()
with open(testdata_fname, 'wb') as fd:
    for chunk in r.iter_content(chunk_size=5096):
        fd.write(chunk)

# unzip dataset
unzip_datasets(download_dir=testdata_dir)

we can either read all bands are define a subset to read

band_selection = ['B02', 'B03', 'B04', 'B05', 'B07', 'B08']

define file-path to ESRI shapefile (all formats understood by fiona work)

import os

base_dir = Path(os.path.dirname(os.path.realpath("file"))).parent.parent in_file_aoi = base_dir.joinpath('data/ref_feats/test_by.shp')

read data from .SAFE dataset for the selected AOI and spectral bands

handler = Sentinel2().from_safe( in_dir=testdata_fname_unzipped, vector_features=in_file_aoi, band_selection=band_selection, apply_scaling=False # if True scales the reflectance values between 0 and 1 ) handler

fig_blue = handler.plot_band('blue', colorbar_label='Surface Reflectance') vals = handler.get_values(['B02'])

m-i-g-u-e-l commented 2 years ago

when switching to numpy<1.22, error disappears.