Reprojection using rioxarray

[Murk89]

I am trying to run a regression with MOD13Q1 NDVI and Weather Research and Forecasting (WRF) climate model outputs.
For this I need to assure that the MODIS raster and WRF dataset are on the same projection and spatial resolution. MODIS geotiff is downloaded in its native projection and resolution using Earth Engine. Then converted the file to a netCDF using gdal package in python.
Next I am trying to use python's reproject_match to reproject MOD13Q1 to WRF's projection, which is explained here: https://spatialreference.org/ref/sr-org/wrf-lambert-conformal-conic/

Following is the code, the current outputs and the desired output.

The code is:

import os
import rioxarray
import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
import rasterio
from rasterio.crs import CRS

import geopandas as gpd
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
from rioxarray.crs import crs_from_user_input
import pandas as pd


#WRF dataset
wrf=xr.open_dataset("wrfout_d02_2007-12-01_00_00_00.nc", decode_coords="all" )

wrf.dims
#returns:
#Dimensions: (Time: 1, south_north: 111, west_east: 114)

#renaming the west_east and south_north dims of wrf to y and x. Because otherwise they are not recognized.
wrf = wrf.rename({'south_north':'y'})
wrf = wrf.rename({'west_east':'x'})

#writing the crs, using proj4
crs = CRS.from_string ('+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
wrf.T2.rio.write_crs("crs", inplace=True)

#MODIS dataset
ndvi=xr.open_dataset("mean_ndvi.nc", decode_coords="all")
ndvi.dims
#returns:
Frozen({'lat': 5949, 'lon': 12720})

#rename lat lon as y and x
ndvi=ndvi.rename({'lon':'x'})
ndvi=ndvi.rename({'lat':'y'})

#crs of the ndvi dataset
ndvi['Band1'].rio.write_crs("epsg:4326", inplace=True)

trial=ndvi.Band1.rio.reproject_match(wrf.T2)
trial.plot()

the above code produces the following plot:

And I actually need it to look like this to be able to do a regression at each x,y

removing the string quotes from "crs" here: wrf.T2.rio.write_crs("crs", inplace=True) generates this plot:

[snowman2]

Are you able to provide the input datasets?

[Murk89]

Yes I can, that's not a problem. But this page doesn't allow for .nc files to upload.

[snowman2]

You should be able to zip it up and then upload.

[Murk89]

@snowman2
I have tried with a compressed zip folder, doesn't allow even then. Can I email the zipped folder to you if you don't mind?

[snowman2]

You may need to upload to another location such as Google drive or Dropbox and share a link.

[Murk89]

The datasets should be accessible here:
https://www.dropbox.com/sh/mmbilrvxo7ly7y7/AAAPYpS7arlSQVYEDz5KJ2y9a?dl=0

[snowman2]

The WRF data has 2D longitude/latitude coordinates.

rioxarray requires 1D x and y variables. Since the coordinates are geographic, you will need to re-project them to the projected coordinate system ('+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs') and flatten it to 1D. See #209.

[snowman2]

Yes, that is correct. However, the x,y data should be coordinates with the dim names as x,y.

[Murk89]

@snowman2
I have added the lat,long transformed x_y to the dataset; have attached the screenshot of the coordinates of the dataset. However the x,y don't come up with any dim names.

Also, then I defined the projection of the new dataset.
Following is the revised code:

import os
import rioxarray
import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
import rasterio
from rasterio.crs import CRS
import geopandas as gpd
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
from rioxarray.crs import crs_from_user_input
import pandas as pd

##opening the WRF dataset:
ds2 = xr.open_dataset("/try.nc")
print (ds2)
##returns:
##<xarray.Dataset>
##Dimensions: (Time: 1, south_north: 111, west_east: 114, x,y: 1)
##Coordinates:
##XLAT (Time, south_north, west_east) float32 ...
##XLONG (Time, south_north, west_east) float32 ...
##XTIME (Time) datetime64[ns] ...
##* x,y (x,y) object 'x_y'
##Dimensions without coordinates: Time, south_north, west_east
##Data variables:
##T2 (x,y, Time, south_north, west_east) float32 ...

crs=CRS.from_string('+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
ds2.T2.rio.write_crs("crs", inplace=True)

##loading the ndvi dataset
ndvi=xr.open_dataset("/mean_ndvi.nc", decode_coords="all")
ndvi=ndvi.rename({'lon':'x'})
ndvi=ndvi.rename({'lat':'y'})

ndvi['Band1'].rio.write_crs("epsg:4326", inplace=True)

##reproject_match
trial=ndvi.Band1.rio.reproject_match(ds2.T2)

##returns: MissingSpatialDimensionError: y dimension not found. 'rio.set_spatial_dims()' or using 'rename()' to change the dimension name to 'y' can address this. Data variable: T2
##end

I think I am missing the step of flattening the x,y of the ds2.T2? I don't understand how to flatten it. Any specific code for that? Many thanks.

[snowman2]

I haven't used the WRF code before. If may be easier to use pyproj.Transformer to transform the 2D lat/lon data and select the first row/column of the coordinates depending on the axis to flatten it.

[Murk89]

Thanks for the all the support @snowman2
I don't understand how I am going to solve this. pyproj is beyond my skills! But thanks anyway. I'll get back if this ever resolves!

[Murk89]

@snowman2
Hi,

I was finally able to reproject the MODIS NDVI onto the WRF grid using salem https://salem.readthedocs.io/en/v0.2.3/installing.html