How can I get Landsat 8 data at southern Latitudes?

[cullen-molitor]

I am trying to use stackstac to get Landsat 8 imagery. This works without issue in any northern latitude, but as soon as I switch to southern latitudes it returns all NA values. After doing a bit of research it appears that Landsat is stored in only northern UTM and any southern zones are given a latitude that is 10,000,000 meters different.

Here is the simplified code I am using to produce this problem. Any help would be greatly appreciated!

import stackstac
import pystac_client
import pyproj
import planetary_computer as pc

# Positive lat values work
# Negative lat values fail
lat =  -5
lon = 20

catalog = pystac_client.Client.open('https://planetarycomputer.microsoft.com/api/stac/v1')

search = catalog.search(
    collections=['landsat-8-c2-l2'],
    intersects=dict(type="Point", coordinates=[lon, lat]),
    datetime="2020-01-01/2020-05-01",
    query={"eo:cloud_cover": {"lt": 10}}
)

items = pc.sign(search)

stack = stackstac.stack(items, assets=["SR_B4", "SR_B3", "SR_B2"])

x_utm, y_utm = pyproj.Proj(stack.crs)(lon, lat)

# I added this based on the link above. This helps subset the image correctly 
# but I am still getting all NA values
if lat < 0:
    y_utm = y_utm - 10000000
else:
    pass

buffer = 2000  # meters

aoi = stack.loc[..., y_utm+buffer:y_utm-buffer, x_utm-buffer:x_utm+buffer]

data = aoi.compute()

I also get a warning after the above chunk

/srv/conda/envs/notebook/lib/python3.8/site-packages/stackstac/prepare.py:413: FutureWarning: pandas.Float64Index is 
deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  xs = pd.Float64Index(np.linspace(minx, maxx, width, endpoint=False))

data.plot.imshow(row="time", rgb="band", robust=True, size=6)

[calebrob6]

Hmmm, stack.crs is epsg:32734 while the data crs (e.g. crs returned by rasterio if you open one of the tiffs) is epsg:32634. If you use x_utm, y_utm = pyproj.Proj("epsg:32634")(lon, lat) then stack.loc behaves as expected (without the need to subtract 10000000).

[cullen-molitor]

Thanks for the fast reply as always! I was hoping to avoid switching CRS manually between location in different zones. Were you able to get an image to plot with this change? I am still getting all NA values and an empty plot.

[calebrob6]

Hmmm, yep -- I get all NA values too.

Two things:

• The data itself isn't messed up
• If I manually set the item to the correct CRS (items[0].properties["proj:epsg"] = 32634) then everything works as expected and I get the following plot
  

@TomAugspurger it looks like the STAC Item here is incorrect?

[TomAugspurger]

Might this be the same issue as stactools-packages/landsat#12? We have an outstanding item to fix those STAC items.

[cullen-molitor]

This does appear to be the same issue. I was able to implement a temporary change based on Caleb's response by using like items[0].properties["proj:epsg"] = items[0].properties["proj:epsg"] - 100 for southern latitudes which takes care of the problem no matter which UTM zone. I am not sure if this is the best solution but it works for now and I will keep an eye on the above issue and hopefully this will be solved soon. I appreciate the help!

[cullen-molitor]

This seems like it is no longer a problem with Landsat. Thanks for the help!

[gilbertocamara]

Hi @cullen-molitor @TomAugspurger, try using the sits R package. It works easily for all hemispheres. An example follows:

# Read a shapefile thar covers the city of Brasilia
shp_file <- system.file("extdata/shapefiles/df_bsb/df_bsb.shp", 
                        package = "sits")
sf_bsb <- sf::read_sf(shp_file)
# select the cube
s2_L8_cube_MSPC <- sits_cube(
        source = "MSPC",
        collection = "LANDSAT-C2-L2",
        roi = sf_bsb,
        platform = "LANDSAT-8", 
        start_date = "2019-06-01",
        end_date = "2019-10-01",
        bands = c("BLUE", "NIR08", "SWIR16", "CLOUD")
)