Slight Difference Between Mc-V and GDAL rectang. projection

[Diego Souza]

Hello all,

I'm trying to visualize GOES-16 imagery (NetCDF) with Mc-V.

When I open a full disk image, it works OK, but it is very slow in my computer.

So I'm trying to subsect the original image using GDAL before opening it with Mc-V.

I have managed to cut the region I wanted with GDAL and open it with Mc-V, however, when overlaying both, there is a slight difference between the original and the cutout.

Please take a look at the attached GIF. The image with the "CMI full disk" legend is the original file. The image with the "gdal test" legend is the cutout.

Please find below the GDAL instructions I'm using to subsect:

To extract the CMI image from the NetCDF and apply a georeference:
gdal_translate -of netCDF -a_srs "+proj=geos +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs +a=6378137. +b=6356752.31414 +lon_0=-89.5 +f+.00335281068119356027489803406172 +h=35786023. +sweep=x" -a_ullr -5433892.8764120 5433892.8764120 5433892.8764120 -5433892.8764120 HDF5:"OR_ABI-L2-CMIPF-M3C09_G16_s20171091800379_e20171091811152_c20171091811228.nc"://CMI GdalTest.nc

To subsect a region and reproject to rectangular:
gdalwarp -of netCDF -te -95 6 -61 24 -s_srs "+proj=geos +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs +a=6378137. +b=6356752.31414 +lon_0=-89.5 +f=.00335281068119356027489803406172 +h=35786023. +sweep=x" -t_srs "+proj=latlong +datum=WGS84" -overwrite GdalTest.nc GdalTest_sub.nc

I know that GDAL is not a tool discussed in the forum, however I'm just adding the code for you to understand what I'm doing.

In the second line of code, the "-s_srs" parameter is the source projection (geostationary) and the "-t_srs" is the target projection (rectangular).

Do you have any idea of what could be generating this difference?

Thank you very much. A salute from Brazil!

[bobc]

Hello -

I don't believe that anyone in our group has enough experience with GDAL to be able to decipher what might be going on. However, you may be able to make use of McIDAS-V's loadGrid() function to quickly display a geographical subset of the data. Here's a script you can run from the Jython Shell (Tools > Formulas > Jython Shell from the Main Display window):

Code: Select all

# specify some variables to point at the ABI data
dataDir = "C:/Users/rcarp/Data/ABI/CMIPF/"
dataFile = "OR_ABI-L2-CMIPF-M3C09_G16_s20171091800379_e20171091811152_c20171091811228.nc"

# dictionary for loadGrid()
abiParms = dict(
    filename = dataDir+dataFile,
    field = "CMI",
    )

# loadGrid call for full domain of data
fullData = loadGrid(**abiParms)

# loadGrid call for partial domain of data
partData = loadGrid(latLonBounds = (24, -95, 6, -61), **abiParms)

# display both layers in the current window
fullLayer = activeDisplay().createLayer("Image Display", fullData)
partLayer = activeDisplay().createLayer("Image Display", partData)

# set enhancement so both layers have same range
fullLayer.setEnhancement("Inverse Gray Scale", range=(170, 270))
partLayer.setEnhancement("Inverse Gray Scale", range=(170, 270))

You will need to change the dataDir and dataFile variables to point towards a file on your machine. To verify that the subsetted and full domain data are the same, I load and display both the full disk as well as the geographical subset of the data. You can toggle the visibility check boxes in the Legend of the Main Display and probe the display with your mouse to verify. If you are only looking to display the subsetted data, you can also remove the fullData and fullLayer lines.

Does this get you what you are looking for without using GDAL?

Thanks -
Bob Carp
McIDAS Help Desk

[Diego_Souza]

Wow! It worked!!! And it is quite fast!

Thanks a lot, I don't need GDAL anymore.

Now I'll study how to apply this to a formula so users do not have to deal with scripts.

Thanks!!!