0

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
1
  • you can also access the projection info with da.spatial_ref.projected_crs_name
    – gman
    Oct 27, 2022 at 16:39

1 Answer 1

2

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

3
  • 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.
    – Rafael
    Oct 27, 2022 at 16:43
  • @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).
    – user2856
    Oct 28, 2022 at 7:49
  • 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
    – Rafael
    Oct 28, 2022 at 14:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.