Add lat and lon to DataArray read in by rioxarray

Question

I downloaded a GeoTIFF from here: https://www.nass.usda.gov/Research_and_Science/Crop_Progress_Gridded_Layers/index.php

(file also available: https://drive.google.com/file/d/1XcfEw-CZgVFE2NJytu4B1yBvjWydF-Tm/view?usp=sharing)

Looking at one of the weeks in 2021, I'd like to add the lat/lon information to the DataArray

I tried:

import rioxarray
fl = 'data/cpc2021/corn/cpccorn2021/condition/cornCond21w24.tif'
da = rioxarray.open_rasterio(fl, masked=True)

However, this returns a DataArray with x and y that don't seem to correspond to the lat/lon. How can I add (or retain) this information?

Expect lat/lon points to be in contiguous US

I found a meta file that has the projections: NAD_1983_Contiguous_USA_Albers

which I believe corresponds to EPSG:5070 (also seen later in the same XML file)

I also found the bounding box for lat/lon coordinates:

 <GeoBndBox esriExtentType="search">
 <exTypeCode Sync="TRUE">1</exTypeCode>
 <westBL Sync="TRUE">-127.360895</westBL>
 <eastBL Sync="TRUE">-68.589171</eastBL>
 <northBL Sync="TRUE">51.723828</northBL>
 <southBL Sync="TRUE">23.297865</southBL>

However, still uncertain how to include this information in my quest to convert to dataframe.

Result of print(da) is:

<xarray.DataArray (band: 1, y: 320, x: 479)>
[153280 values with dtype=float32]
Coordinates:
 * band (band) int64 1
 * x (x) float64 -2.305e+06 -2.296e+06 ... 1.987e+06 1.996e+06
 * y (y) float64 3.181e+06 3.172e+06 ... 3.192e+05 3.102e+05
 spatial_ref int64 0
Attributes:
 AREA_OR_POINT: Area
 RepresentationType: ATHEMATIC
 STATISTICS_COVARIANCES: 0.1263692188822515
 STATISTICS_MAXIMUM: 4.8569073677063
 STATISTICS_MEAN: 3.7031858480518
 STATISTICS_MINIMUM: 2.1672348976135
 STATISTICS_SKIPFACTORX: 1
 STATISTICS_SKIPFACTORY: 1
 STATISTICS_STDDEV: 0.35548448472789
 scale_factor: 1.0
 add_offset: 0.0

• gman: you can also access the projection info with da.spatial_ref.projected_crs_name

Answer 1

Your data is provided in the NAD_1983_Contiguous_USA_Albers projection (just check gdalinfo <filename.tif> in the console). If you want it in Longitude/Latitude, you'll need to reproject it to latitude/longitude (the target spatial reference, given by EPSG code 4326. This will work on a single file:

import rioxarray
_da = rioxarray.open_rasterio("cornProg22w30.tif", masked=True)
da = da.rio.reproject("EPSG:4326")

This gives

<xarray.DataArray (band: 1, y: 268, x: 500)>
array([[[nan, nan, nan, ..., nan, nan, nan],
 [nan, nan, nan, ..., nan, nan, nan],
 [nan, nan, nan, ..., nan, nan, nan],
 ...,
 [nan, nan, nan, ..., nan, nan, nan],
 [nan, nan, nan, ..., nan, nan, nan],
 [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
Coordinates:
 * x (x) float64 -117.3 -117.2 -117.1 ... -69.87 -69.77 -69.68
 * y (y) float64 49.45 49.35 49.26 49.16 ... 24.24 24.14 24.04 23.95
 * band (band) int64 1
 spatial_ref int64 0
Attributes:
 RepresentationType: ATHEMATIC
 STATISTICS_COVARIANCES: 0.01716971322246022
 STATISTICS_MAXIMUM: 0.87459498643875
 STATISTICS_MEAN: 0.48183899334969
 STATISTICS_MINIMUM: 0.32959023118019
 STATISTICS_SKIPFACTORX: 1
 STATISTICS_SKIPFACTORY: 1
 STATISTICS_STDDEV: 0.13103325235397
 scale_factor: 1.0
 add_offset: 0.0

and looks like this enter image description here

• Rafael: thanks! where did you find that NAD_1983_Contiguous_USA_Albers corresponds to EPSG:4326? In the XML file that came with the data as well as some quick searches I thought it would be EPSG:5070.
• user2856: @Rafael NAD_1983_Contiguous_USA_Albers is not EPSG:4326 it is EPSG:5070 which uses metres as the horizontal unit (x & y xoordinates). The code in the answer reprojects the data from EPSG:5070 to EPSG:4326 which uses decimal degrees (lat & lon coordinates).
• Rafael: I asked something similar here: stackoverflow.com/questions/74208825/… will post your solution as an answer but if you'd like to post there as well I'll happily accept your answer on stackoverlow as well