[Bug]: Getting incorrect number of pixels and their X/Y coordinates or incorrect data when working in non-EPSG:4326 coordinates

[shahriar49]

Bug Summary

I want to get an area from a Landsat Collection under Google Earth Engine. I set the coordinates of bounding box of 333x333 pixels (the exact numbers that match the UTM pixel bounds) one time in the UTM Zone 13 values and another time in corresponding lat/lon values and call the xr.open_dataset() function with projection set to ee.Projection(crs='EPSG:32613', transform=[30, 0, 699960, 0, -30, 4300020]). In both cases it throws an error

Traceback (most recent call last):
  File "C:\Python3.9.13\lib\contextlib.py", line 137, in __exit__
    self.gen.throw(typ, value, traceback)
  File "C:\Users\sheydari\Python39-2\lib\site-packages\xarray\core\utils.py", line 628, in close_on_error
    yield
  File "C:\Users\sheydari\Python39-2\lib\site-packages\xee\ext.py", line 1153, in open_dataset
    ds = store_entrypoint.open_dataset(
  File "C:\Users\sheydari\Python39-2\lib\site-packages\xarray\backends\store.py", line 43, in open_dataset
    vars, attrs = filename_or_obj.load()
  File "C:\Users\sheydari\Python39-2\lib\site-packages\xarray\backends\common.py", line 221, in load
    (_decode_variable_name(k), v) for k, v in self.get_variables().items()
  File "C:\Users\sheydari\Python39-2\lib\site-packages\xee\ext.py", line 676, in get_variables
    vars_ = [(name, self.open_store_variable(name)) for name in self._bands()]
  File "C:\Users\sheydari\Python39-2\lib\site-packages\xee\ext.py", line 676, in <listcomp>
    vars_ = [(name, self.open_store_variable(name)) for name in self._bands()]
  File "C:\Users\sheydari\Python39-2\lib\site-packages\xee\ext.py", line 566, in open_store_variable
    arr = EarthEngineBackendArray(name, self)
  File "C:\Users\sheydari\Python39-2\lib\site-packages\xee\ext.py", line 791, in __init__
    x_size = max(1, int(np.round((x_max - x_min) / np.abs(self.store.scale_x))))
ValueError: cannot convert float NaN to integer

If I set the projection argument to the projection object that I get from the Landsat image itself, it works but returns a dataset with incorrect number of rows and columns (several more rows and one less column) and incorrect pixel center locations (they are all off by several meters in both directions).

If I give the xr.open_dataset() function the crs and scale separately (not through the projection argument) it works in setting the dataset dimension and coordination values correctly. However, it returns the dimension names as lat and lon and crs as EPSG:4326 in xarray dataset, which is incorrect, and the dataset variables (Landsat band values) are all null.

P.S. It also works only with 'EPSG:32613' setting for crs, not with the wkt string I get from GEE's Landsat projection.

I appreciate if you look at this issue promptly.

Steps to Reproduce

You can use this code to see the behavior.

point = ee.Geometry.Point([-102.05, 38.05])
image = ee.ImageCollection("NASA/HLS/HLSL30/v002").filterDate("2018-08-01", "2018-09-01").filterBounds(point).select(['B4']).first()
projection = image.projection()
crs = 'EPSG:32613'
crsTransform = projection.getInfo()['transform']
scale = projection.nominalScale().getInfo()
BBox1 = ee.Geometry.Polygon([[[750000, 4210020],[759990, 4210020],[759990, 4220010],[750000, 4220010]]], crs, evenOdd=False)
BBox2 = ee.Geometry.Polygon([[[-102.15268713055727, 38.0034113624073], [-102.03904576692914, 38.0006026911579],
    [-102.0354189991692, 38.09052020120823],[-102.14919908937703, 38.09333791178974]]], 'EPSG:4326', evenOdd=False)
ds = xr.open_dataset(
    ee.ImageCollection([image]),
    engine='ee',
    geometry=BBox1,
    # projection=ee.Projection(crs=crs,transform=crsTransform)
    projection=projection,
    # projection=crs,
    # scale=scale
)
print(ds)

Current Behavior

Explained above.

Expected Behavior

Get the correct dimension, coordination names, and dataset values.

Relevant log output

Xee Version

0.0.19

Contact Details

No response

Code of Conduct

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

[DecaiJin]

@shahriar49 I get the same error when the xr.open_dataset accept the projection. After checking source code, I found that ee and pyproj are used in combination, resulting in two reprojection of "geometry" (#207 and #217 already mentions this issue). In addtion, if you only provide "crs" to open_dataset, the reprojection of pyproj is not accurate compared to ee, which is the ultimate cause of dimensional mismatch and center coordinate displacement. By the way, I hope the authors can rewrite open_dataset function #217

[DecaiJin]

@shahriar49 The temporary solution is to comment out xee/ext.py line 676-677

[shahriar49]

Thank you very much @DecaiJin for your attention. But the line 676 in my xee/ext.py is:

vars_ = [(name, self.open_store_variable(name)) for name in self._bands()]

and line 677 is just a blank line. vars_ variable is used in line 679.

Also, in the code that is currently on the git, the line 676 is:

return x_min, y_min, x_max, y_max

Do you mean comment out this? Am I looking at the wrong file?

[DecaiJin]

@shahriar49 Sorry, may be I have modified the code before check lines. The right line is Line 675-675 on the git. It can be replaced by:
x_min, y_min = (x_min_0, y_min_0)
x_max, y_max = (x_max_0, y_max_0)
You can provide the original projection to open_dataset. The "crs" doesn't matter anymore.

[shahriar49]

Thank you very much @DecaiJin for your help. I changed the lines and the result depends on what I pass to the xr.open_dataset function as projection:

1- If I set projection=ee.Projection(crs='EPSG:32613',transform=crsTransform), the export is done correctly. I mean, the exported projection is EPSG:32613 and dataset has X values at (750015.0, 750045.0, ..., 759975.0) and Y values at (4210035.0, 4210065.0, ..., 4219995.0) with shapes of 333, 333. The original bound was at X:(750000 ~ 759990) and Y:(4210020 ~ 4220010).

2- If I set projection=image.projection(), the export is wrong. It returns just one pixel at a totally wrong lat/lon coordinate with projection of EPSG:4326.

3- If I set projection='EPSG:32613', scale=30, exported X and Y values are correct but the projection is EPSG:4326 and all of the exported values are NaN.

4- If I I set projection=ee.Projection(crs=image.projection().getInfo()['wkt'],transform=crsTransform), it returns just one pixel at a totally wrong lat/lon coordinate but claims that the projection is UTM Zone 13.

5- If I set projection=image.projection().getInfo()['wkt'], if returns just one pixel somewhat close to the area boundary with projection of EPSG:4326.

So providing the original projection (case 2 above) does not work, which is very strange to me. I can go ahead with option 1, but the catch is that I should give it the crs as a well-known EPSG code and neither the ee object nor the wkt string do not work.

If I just want to continue with the fix you mentioned and use the first of five approaches I mentioned above, will there be any other issue that can possibly rise somewhere else in the code?

[DecaiJin]

@shahriar49 I did not figure out why you are having the problem with scenario 2. I do not have that problem in my example. However, I'm different from you: projection=image.select("band_name").projection()

[shahriar49]

@DecaiJin Thanks again for reply. I have already selected my band before doing calling the open_dataset. I tried my code again and my scenario 2 above still returns a 1 pixel x 1 pixel output.

[Lkruitwagen]

@shahriar49 I think you were really close. I think the geometry needs to be passed in as a geographic geography (i.e. latlon/epsg:4326). And the the projection needs to have your target (i.e. UTM) projection, as well as the transform, which you image here has.

This worked fine for me:

point = ee.Geometry.Point([-102.05, 38.05])
    image = (
        ee.ImageCollection("NASA/HLS/HLSL30/v002")
        .filterDate("2018-08-01", "2018-09-01")
        .filterBounds(point)
        .select(['B4'])
        .first()
    )
    projection = image.projection()

    ee_bbox_4326 = ee.Geometry.Polygon([[
        [-102.15268713055727, 38.0034113624073], 
        [-102.03904576692914, 38.0006026911579],
        [-102.0354189991692, 38.09052020120823],
        [-102.14919908937703, 38.09333791178974]
    ]], 'EPSG:4326', evenOdd=False)
    ds = xr.open_dataset(
        ee.ImageCollection([image]),
        engine='ee',
        geometry=ee_bbox_4326,
        projection=projection,
    )

[Lkruitwagen]

Or rather, squeeze and transpose:

[shahriar49]

Thanks @Lkruitwagen for your attention but the result is not right. The result should be a 333x333 array, not a 332x354 array, and the pixel centers are not at the right locations (15 meters after the box corner). I checked again and see that the easting/northing values that I provided (which show a 333x333 @ 30m resolution area) exactly map to the above lat/lon values. So the bounding box in lat/lon is correct but the xee module is not doing the job right.
Previously I was using xee version 0.0.19 and I was able to get the correct dataset by passing the bounding box in easting/northing values and passing the projection as I specified in my above post dated Mar 31. I see the current version is 0.0.20. So I upgraded my xee to that version but the problem persists.
I appreciate if you - or other colleagues - can help resolve the issue. Without the guarantee for proper conversion, this module will be useless.