[Bug]: EEException: Unable to transform edge

[DavidCarMor]

Bug Summary

I am creating a day and night image collection from the ERA5 hourly image collection, and it appears that xarray works until I attempt to plot one of the images to verify that the code is functioning correctly. It looks like there is an error on the borders of the images, despite I clipped the images to avoid that issue. Has anybody faced a similar issue?

Steps to Reproduce

Here is an example of the code that I am using.

import ee
import math
import xarray
ee.Initialize(opt_url='https://earthengine-highvolume.googleapis.com')

# Function to compute sunrise/sunset hours from image date and coordinates
def getSunrise(image):
    PI = ee.Number(math.pi)
    # Use Earth Engine's pixelLonLat to get latitude and longitude bands
    coords = ee.Image.pixelLonLat()
    lat = coords.select('latitude')
    lon = coords.select('longitude')
    # Extract date information
    date = image.date()
    year = ee.Number.parse(date.format('y'))
    doy = ee.Number.parse(date.format('D'))
    # Total days in the year
    days_in_year = ee.Date.fromYMD(year.add(1), 1, 1).difference(ee.Date.fromYMD(year, 1, 1), 'day')
    # Fractional year in radians (gamma)
    gamma = doy.subtract(1).multiply(2 * math.pi).divide(days_in_year)
    # Equation of time (in minutes)
    eq_time = ee.Number(229.18).multiply(
        ee.Number(0.000075)
        .add(ee.Number(0.001868).multiply(gamma.cos()))
        .subtract(ee.Number(0.032077).multiply(gamma.sin()))
        .subtract(ee.Number(0.014615).multiply(gamma.multiply(2).cos()))
        .subtract(ee.Number(0.040849).multiply(gamma.multiply(2).sin()))
    )
    # Solar declination (in radians)
    decl = ee.Number(0.006918) \
        .subtract(ee.Number(0.399912).multiply(gamma.cos())) \
        .add(ee.Number(0.070257).multiply(gamma.sin())) \
        .subtract(ee.Number(0.006758).multiply(gamma.multiply(2).cos())) \
        .add(ee.Number(0.000907).multiply(gamma.multiply(2).sin())) \
        .subtract(ee.Number(0.002697).multiply(gamma.multiply(3).cos())) \
        .add(ee.Number(0.00148).multiply(gamma.multiply(3).sin()))
    # Latitude in radians
    lat_rad = lat.multiply(PI.divide(180))
    # Hour angle (in radians)
    cos_ha = ee.Image.constant(math.cos(math.radians(90.833))) \
        .divide(lat_rad.cos().multiply(decl.cos())) \
        .subtract(lat_rad.tan().multiply(decl.tan()))
    # Clip cosH to [-1, 1] to avoid NaNs
    ha = cos_ha.clamp(-1, 1).acos()
    # Compute sunrise and sunset times (in minutes from midnight UTC)
    sunrise_min = ee.Image.constant(720) \
        .subtract(lon.add(ha.multiply(ee.Number(180).divide(PI))).multiply(4)) \
        .subtract(eq_time)
    
    sunset_min = ee.Image.constant(720) \
        .subtract(lon.subtract(ha.multiply(ee.Number(180).divide(PI))).multiply(4)) \
        .subtract(eq_time)
        
    # Convert to hours and cast to Uint8
    sunrise_hour = sunrise_min.divide(60).toUint8().rename('sunriseHour')
    sunset_hour = sunset_min.divide(60).toUint8().rename('sunsetHour')

    return sunrise_hour.addBands(sunset_hour)


def getRH(year):
  start = ee.Date.fromYMD(year, 1, 1) #inclusive date
  end = start.advance(1, 'year') # exclusive date
  ERA5 = ee.ImageCollection("ECMWF/ERA5_LAND/HOURLY").filterDate(start,end)\
    .select(['dewpoint_temperature_2m','temperature_2m','surface_pressure'])
  
  # Getting ERA5 projection
  ERA5projection = ERA5.first().projection()
  ERA5crs = ERA5projection.crs().getInfo()
  ERA5scale = ERA5projection.nominalScale().getInfo()
  
  # Get sunrise/sunset image from first hourly image
  sunTimes = getSunrise(ERA5.first()).reproject(ERA5crs,None,ERA5scale)
  
  # List of Days of the Year (from 0 to 364 or 365 for leap years)
  DoY = ee.List.sequence(0, end.difference(start, 'days').subtract(1))
  
  def map_day(day_offset):
    date = start.advance(day_offset, 'day')
    DP = ERA5.filterDate(date, date.advance(1, 'day'))

    def add_hour(img):
        hour = ee.Number.parse(img.date().format('H')).toUint8()
        return img.addBands(ee.Image.constant(hour).rename('hour').toUint8())

    DP = DP.map(add_hour)

    def mask_day(img):
        hour = img.select('hour')
        mask = hour.gte(sunTimes.select('sunriseHour')).And(hour.lte(sunTimes.select('sunsetHour')))
        return img#.updateMask(mask.Or(img.mask()))

    def mask_night(img):
        hour = img.select('hour')
        mask = hour.lt(sunTimes.select('sunriseHour')).Or(hour.gt(sunTimes.select('sunsetHour')))
        return img#.updateMask(mask.Or(img.mask()))

    day_imgs = DP.map(mask_day)
    night_imgs = DP.map(mask_night)

    date_img = (day_imgs.select(['dewpoint_temperature_2m', 'temperature_2m', 'surface_pressure'])
                    .mean().rename(['Td_day', 'Ta_day', 'P_day'])
                    .setDefaultProjection(ERA5crs, None, ERA5scale)
                    .addBands(night_imgs.select(['dewpoint_temperature_2m', 'temperature_2m', 'surface_pressure'])
                              .mean().rename(['Td_night', 'Ta_night', 'P_night'])
                              .setDefaultProjection(ERA5crs, None, ERA5scale))
                    .set('system:time_start', date.millis()))

    # ERA5 air and dew point temperatures at 2 meters
    Ta = date_img.select(['Ta_day', 'Ta_night'])
    Td = date_img.select(['Td_day', 'Td_night'])
    TaC = Ta.subtract(273.15)
    TdC = Td.subtract(273.15)

    # Masks
    ice_mask = Ta.lte(ee.Image.constant(250.16)).rename(['ice_day', 'ice_night'])
    mix_mask = Ta.gt(ee.Image.constant(250.16)).And(Ta.lt(ee.Image.constant(273.16))).rename(['mix_day', 'mix_night'])

    # Teten's (ECMWF) water vapor pressure at saturation over water and ice [Pa]
    # esat with respect to ice (i) and water (w)
    easat_w = Ta.subtract(273.16).divide(Ta.subtract(32.19)).multiply(17.502).exp().multiply(611.21)
    easat_i = Ta.subtract(273.16).divide(Ta.subtract(-0.7)).multiply(22.587).exp().multiply(611.21)

    # mixed phase parameter
    alpha = Ta.subtract(250.16).divide(23).pow(2)\
      .where(ice_mask,0).where(mix_mask.Or(ice_mask).Not(),1)

    # Teten's (ECMWF) water vapor pressure (weighted average) [Pa]
    easat = alpha.multiply(easat_w).add(ee.Image.constant(1).subtract(alpha).multiply(easat_i)).rename(['easat_day','easat_night'])

    # pv
    pv = Td.multiply(17.67).exp().multiply(611.21).divide(Td.add(243.5)).rename(['e_day', 'e_night'])

    # RHte
    rh_tt = pv.divide(easat).clamp(0, 1).rename(['RHtt_day', 'RHtt_night'])

    # RHg
    rh_g = (TdC.multiply(17.625).divide(Td.add(243.04-273.16)).exp()
            .divide(TaC.multiply(17.625).divide(Ta.add(243.04-273.16)).exp())
            .clamp(0, 1).rename(['RHg_day', 'RHg_night']))

    return date_img.addBands([rh_tt,rh_g]).set('system:time_start', date.millis())\
                              .reproject(ERA5crs, None, ERA5scale)

  yearDP = ee.ImageCollection(DoY.map(map_day))

  return yearDP
  
  startYear = 2012
endYear = 2013
yearList = ee.List.sequence(startYear, endYear).getInfo()

lclIC = getRH(startYear)

# Define a bounding box from -60 to +60 degrees latitude and global longitude
region = ee.Geometry.Rectangle([-180, -90, 180, 90])

# Apply the clipping to each image in the ImageCollection
lclIC = lclIC.map(lambda img: img.clip(region))

IC_proj = lclIC.first().select('Ta_day').projection().getInfo()
IC_crs = IC_proj['crs']
IC_scale = IC_proj['transform'][0]
ds = xarray.open_dataset(lclIC.filterDate('2012-01-01','2012-01-03'), engine='ee', crs=IC_crs, scale=IC_scale)

data_test = ds.isel(time=0).copy()

import matplotlib.pyplot as plt
import numpy as np

# Plot the LCL variable
plt.figure(figsize=(10, 6))
plt.pcolormesh(
    data_test["Ta_day"].lon,
    data_test["Ta_day"].lat,
    data_test["Ta_day"].T, # 'LCLrte_night', 'RHg_night', RHte_night, P_night, Ta_night
    cmap="viridis",
    shading="auto"
)
plt.title('Example LCL [m]'), plt.xlabel("Longitude"), plt.ylabel("Latitude")

# Add a colorbar
cbar = plt.colorbar(label="Lifting Condensation Level (meters)")  # Add a colorbar for clarity

# Set tick locations and labels
# cbar.set_ticks(range(0,60,10))

plt.show()

Current Behavior

I am unable to use the images at their native scale. I can only use them if I let Xee reproject them into a 1-degree scale.

Expected Behavior

I should be able to use the images at their native scale.

Relevant log output

EEException                               Traceback (most recent call last)
Cell In[15], [line 6](vscode-notebook-cell:?execution_count=15&line=6)
      [4](vscode-notebook-cell:?execution_count=15&line=4) # Plot the LCL variable
      [5](vscode-notebook-cell:?execution_count=15&line=5) plt.figure(figsize=(10, 6))
----> [6](vscode-notebook-cell:?execution_count=15&line=6) plt.pcolormesh(
...
    [406](https://vscode-remote+olympus-002edeac-002ewfu-002eedu.vscode-resource.vscode-cdn.net/deac/egr/lowmanGrp/miniconda3/envs/gee/lib/python3.11/site-packages/ee/data.py:406)   return call.execute(num_retries=num_retries)
    [407](https://vscode-remote+olympus-002edeac-002ewfu-002eedu.vscode-resource.vscode-cdn.net/deac/egr/lowmanGrp/miniconda3/envs/gee/lib/python3.11/site-packages/ee/data.py:407) except googleapiclient.errors.HttpError as e:
--> [408](https://vscode-remote+olympus-002edeac-002ewfu-002eedu.vscode-resource.vscode-cdn.net/deac/egr/lowmanGrp/miniconda3/envs/gee/lib/python3.11/site-packages/ee/data.py:408)   raise _translate_cloud_exception(e)

EEException: Unable to transform edge (-1799.999861, 644.000000 to -1799.999861, 644.000000) from EPSG:4326 PLANAR [0.10000000000000002, 0.0, 0.0, 0.0, -0.10000000000000002, 0.0] to EPSG:4326.

Xee Version

0.0.20

Contact Details

carcpd21@wfu.edu

Code of Conduct

• I agree to follow this project's Code of Conduct