Search code examples
pythonprojectionpython-xarraypyproj

reproject IMERG using rioxarray

I would like to make reprojected array using rioxarray.

My data is IMERG from NASA and GK-2A weather satellite image. I want to convert Imerg data to a projection of an image in gk-2a.

IMRG files can be obtained from NASA, but an account is required. https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGHH.06/2021/001/3B-HHR.MS.MRG.3IMERG.20210101-S000000-E002959.0000.V06B.HDF5

GK-2A data link: https://nmsc.kma.go.kr/homepage/html/satellite/viewer/selectImgDown.do?fileKey=GK2A:AMI:LE1B:IR105:NC:EA:020:LC&observationTime=202108080530&type=NC

This file is epsg:4326 and I want to change projection to Lambert conformal.

Parameters: 
projection units:   meters       ==> +units=m
datum (spheroid):   WGS_84       ==> +datum=WGS84
1st standard parallel:  30 deg N ==> +lat_1=30
2nd standard parallel:  60 deg N ==> +lat_2=60
Central meridian:   126 deg E    ==> +lon_0=126
Latitude of origin: 38 deg N     ==> +lat_0=38

Area of Interest

upper left
Easting: -2999000m/76.811834E
Northing: 2599000m/53.303712N

upper right
Easting: 2999000m/175.188166E
Northing: 2599000m/53.303712N

Bottom left
Easting: -2999000m/101.395259E
Northing: -2599000m/11.308528N

Bottom right
Easting: 2999000m/150.604741E
Northing: -2599000m/11.308528N

This is code what I tried. (Using xarray interp to reproject a dataarray?)

import xarray as xr
import numpy as np
from pyproj import CRS
import rioxarray
import netCDF4
import matplotlib.pyplot as plt
file_path = 'E:/IMERG/3B-HHR.MS.MRG.3IMERG.20210827-S000000-E002959.0000.V06B.HDF5'



ncf = netCDF4.Dataset(file_path, diskless=True, persist=False)
nch = ncf.groups.get('Grid')
xds = xr.open_dataset(xr.backends.NetCDF4DataStore(nch))
ds = xr.decode_cf(xds,decode_coords="all")
ds = ds.rio.write_crs("EPSG:4326", inplace=True)
ds = ds.precipitationCal.squeeze('time')
ds = ds.transpose('lat', 'lon')

sat_path = 'D:/EA/IR105/20210809/gk2a_ami_le1b_ir105_ea020lc_202108090000.nc'
crs = '+proj=lcc +datum=WGS84 +lat_1=30 +lat_2=60 +lon_0=126 +lat_0=38 +bounds=(-2599000,2599000,-2999000,2999000)'
xds = xr.open_dataset(sat_path)
sat = xds.rio.write_crs(crs, inplace=True)
xds = xds.squeeze().rename_dims({"dim_x": "x", "dim_y": "y"}).transpose('y', 'x')
xds.coords["x"] = xds.x
xds.coords["y"] = xds.y


xds_repr_match = ds.rio.reproject_match(xds)
print(ds)
print(xds_repr_match)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
xds_repr_match.plot(ax=axes[0])
ds.plot(ax=axes[1])

plt.show()

what I got image I got this image but I expect below image(I draw with cartopy).

expected image

Solution

  • This is the part you are missing:

    sat_path = "gk2a_ami_le1b_ir105_ea020lc_202108080000.nc"
crs = "+proj=lcc +datum=WGS84 +lat_1=30 +lat_2=60 +lon_0=126 +lat_0=38"
# no +bounds=(-2999000,2999000,-2599000,2599000) parameter
# https://proj.org/operations/projections/lcc.html

sat_xds = xarray.open_dataset(
        sat_path,
        decode_coords="all",
).image_pixel_values.squeeze().rename_dims(
    {"dim_x": "x", "dim_y": "y"}
).transpose('y', 'x')
sat_xds.rio.write_crs(crs, inplace=True)

    At this point, there are no x or y coordinates on your dataset, so you need to figure out the transform:

    # https://rasterio.readthedocs.io/en/latest/api/rasterio.transform.html#rasterio.transform.from_bounds
transform = rasterio.transform.from_bounds(
    west=sat_xds.upper_left_easting,
    south=sat_xds.lower_right_northing,
    east=sat_xds.lower_right_easting,
    north=sat_xds.upper_left_northing,
    width=sat_xds.image_width,
    height=sat_xds.image_height,
)
sat_xds.rio.write_transform(transform, inplace=True)

    After that, this appears to work:

    xds_repr_match = ds.rio.reproject_match(sat_xds)

    image image