Как сделать массив временных рядов Numpy Изображение в Google Earth Engine Python API - PullRequest
Новая
       9

Как сделать массив временных рядов Numpy Изображение в Google Earth Engine Python API

0 голосов
/ 02 марта 2020

Я работаю с Google EarthEngine (Python API, но это не имеет большого значения). У меня есть коллекция ImageCollection, которую нужно уменьшить на регион для каждого отдельного изображения (изображения массива np) в коллекции.

Я хочу построить коллекцию изображений временного ряда в виде изображения массива np (например, файла NetCDF) для дальнейший анализ. Я построил один код здесь, который я использую для l oop для создания коллекции изображений массива временных рядов np (широта, долгота, время, данные), но я не понял, как это сделать с соответствующим кодом.

Показывает проблему измерения нежелательной ошибки.

Пожалуйста, помогите разработчику решить эту проблему или поищите другие способы / хитрости, а затем предложите мне. 

# -*- coding: utf-8 -*-
"""
Created on Mon Mar  2 15:49:17 2020

@author: HONEY
"""

##------Import Library---------------------------------------------------------
## processing -->>

import time
from datetime import datetime as dt
import warnings
warnings.filterwarnings('ignore')
np.set_printoptions(suppress=True)

%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib
matplotlib.rc('font', **{'size': 16})

# Import the Earth Engine Python Packages
import ee
import ee.mapclient 
ee.Initialize()
##-----------------------------------------------------------------------------



# Define the area
area = ee.Geometry.Polygon([[[-81.76903193287582,27.802143094663233],\
                             [-80.7129711418602,27.802143094663233],\
                             [-80.7129711418602,28.551464360900866],\
                             [-81.76903193287582,28.551464360900866],\
                             [-81.76903193287582,27.802143094663233]]])



Bands = ['BRDF_Albedo_Parameters_vis_iso', 'BRDF_Albedo_Parameters_vis_vol',\
         'BRDF_Albedo_Parameters_vis_geo']

collection1 = ee.ImageCollection('MODIS/006/MCD43A1')\
                      .filterBounds(area)\
                      .filterDate('2018-05-01','2018-05-3')\
                      .select(Bands)


print(" number of image: ",collection1.size().getInfo())

# Get list of images which correspond with the above
images = [item.get('id') for item in collection1.getInfo().get('features')]


# Loop over all images and make collection of np array image
for image in images:
    print('Albedo_Data: ',image)
    im = ee.Image(image)
#    CRS=im.projection().getInfo()['crs']
    CRS='EPSG:4326'
    # Obtain date from timestamp in metadata
    date = dt.fromtimestamp(im.get("system:time_start").getInfo() / 1000.)
    print(date,"date")
    # Extract pixel value
    data= im.select(Bands)
#   print(data.getInfo())
    latlon = ee.Image.pixelLonLat().addBands(data)
    countRes=latlon.reduceRegion(reducer=ee.Reducer.count(),\
                        crs=CRS,geometry=collection.geometry().bounds(),\
                        maxPixels=1e50, scale=500)
    dictt = latlon.reduceRegion(ee.Reducer.toList(),\
                                crs=CRS,geometry=collection.geometry().bounds(),\
                                maxPixels=1e50, scale=500)

    lats = np.array((ee.Array(dictt.get("latitude")).getInfo()))
    print("Lats_value:",lats,"Lats_Size:",lats.size)
    lons = np.array((ee.Array(dictt.get("longitude")).getInfo()))
    print("lons_value:",lats,"lons_Size:",lats.size)
    #print("First five (lon,lat) pairs:"); list(zip(lons[:5], lats[:5]))
    #variables(Bands)
    dim=np.sqrt(countRes.get('BRDF_Albedo_Parameters_vis_iso').getInfo()).astype(np.int16)
    vis_iso = np.array((ee.Array(dictt.get("BRDF_Albedo_Parameters_vis_iso")).getInfo())).astype(np.float)
    vis_vol = np.array((ee.Array(dictt.get("BRDF_Albedo_Parameters_vis_vol")).getInfo())).astype(np.float)
    vis_geo = np.array((ee.Array(dictt.get("BRDF_Albedo_Parameters_vis_geo")).getInfo())).astype(np.float)

    nrows, ncols = dim, dim 


    # plot
    fig = plt.figure(figsize=(16,5))

    ax1 = fig.add_subplot(1,4,1)
    ax1.imshow(lon2d, cmap=cm.plasma); 
    ax1.set_title("x coordinate")

    ax2 = fig.add_subplot(1,4,2, sharey=ax1)
    ax2.imshow(lon2d, cmap=cm.cividis); 
    ax2.set_title("longitude")
...