Using GDALwarp for reprojecting netCDF file?

Question

I am trying to re-project a netCDF file using the gdalwarp command. However, the gdalwarp is giving me a black output.

Below is my code:

import gdal
import osr
import netCDF4
from osgeo.gdalconst import *
import os
import subprocess

gdal.AllRegister()
file_nc = 'input.nc'
ds = gdal.Open(file_nc, GA_ReadOnly)
Dataset = ds.GetSubDatasets()

substr = 'Rrs_655'
band = ''
names = [i[0] for i in Dataset]
for word in names[:]:
if word.endswith(substr):
band = word
b = 'gdal_translate -a_scale 0.000002 -a_offset 0.05 -a_nodata -32767 '+ band + ' output13.tif'
print(b)
subprocess.run(b, shell = True)
cmd = 'gdalwarp output13.tif outfile2.tif -t_srs "+proj=longlat +ellps=WGS84"'
subprocess.run(cmd, shell = True)

The gdalinfo of output13.tif band is:

Files: output13.tif
Size is 7821, 7931
Coordinate System is `'
GCP Projection =
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 7931.0)
Upper Right ( 7821.0, 0.0)
Lower Right ( 7821.0, 7931.0)
Center ( 3910.5, 3965.5)
Band 1 Block=7821x1 Type=Int16, ColorInterp=Gray
NoData Value=-32767
Offset: 0.05, Scale:2e-06
Metadata:
geophysical_data_Rrs_655_add_offset=0.050000001
geophysical_data_Rrs_655_long_name=Remote sensing reflectance at 655 nm
geophysical_data_Rrs_655_scale_factor=2e-06
geophysical_data_Rrs_655_solar_irradiance=1550.3799
geophysical_data_Rrs_655_standard_name=surface_ratio_of_upwelling_radiance_emerging_from_sea_water_to_downwelling_radiative_flux_in_air
geophysical_data_Rrs_655_units=sr^-1
geophysical_data_Rrs_655_valid_max=25000
geophysical_data_Rrs_655_valid_min=-30000
geophysical_data_Rrs_655__FillValue=-32767

The gdalinfo if the file outfile2.tif is:

Files: outfile2.tif
Size is 12259, 12307
Coordinate System is:
GEOGCS["WGS 84",
DATUM["unknown",
SPHEROID["WGS84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]]
Origin = (-20696.575935146218399,3354.169516545325678)
Pixel Size = (1.950162663734099,-1.950162663734099)
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -20696.576, 3354.170) (Invalid angle,Invalid angle)
Lower Left ( -20696.576, -20646.482) (Invalid angle,Invalid angle)
Upper Right ( 3210.468, 3354.170) (Invalid angle,Invalid angle)
Lower Right ( 3210.468, -20646.482) (Invalid angle,Invalid angle)
Center ( -8743.054, -8646.156) (Invalid angle,Invalid angle)
Band 1 Block=12259x1 Type=Int16, ColorInterp=Gray
NoData Value=-32767
Metadata:
geophysical_data_Rrs_655_add_offset=0.050000001
geophysical_data_Rrs_655_long_name=Remote sensing reflectance at 655 nm
geophysical_data_Rrs_655_scale_factor=2e-06
geophysical_data_Rrs_655_solar_irradiance=1550.3799
geophysical_data_Rrs_655_standard_name=surface_ratio_of_upwelling_radiance_emerging_from_sea_water_to_downwelling_radiative_flux_in_air
geophysical_data_Rrs_655_units=sr^-1
geophysical_data_Rrs_655_valid_max=25000
geophysical_data_Rrs_655_valid_min=-30000
geophysical_data_Rrs_655__FillValue=-32767

By the corner coordinates ( Upper Left ( 0.0, 0.0) and so on) output13.tif it is not really georeferenced. That should be fixed first. — user30184, Apr 17 '19 at 18:45
I am not so familiar with netCDF and without test data I can't do anything. Try with gdalinfo and see if even it can read the georeferencing from your subdatasets. — user30184, Apr 17 '19 at 19:13
Here is the link to the input.nc file: https://drive.google.com/file/d/1xw6OB0oZhOde0Lb4x3YJLe4hyyDfO2vu/view?usp=sharing
The commands that I am using are:

gdal_translate -a_scale 0.000002 -a_offset 0.05 -a_nodata -32767 HDF5:"input.nc"://geophysical_data/Rrs_655 output13.tif

gdalwarp output13.tif outfile2.tif -t_srs "+proj=longlat +ellps=WGS84"

Thank you! — npal, Apr 17 '19 at 19:25
You seem to have HDF5 data and driver page https://www.gdal.org/frmt_hdf5.html says There is no universal way of storing georeferencing in HDF files. However, some product types have mechanisms for saving georeferencing, and some of these are supported by GDAL. Currently supported are (subdataset_type shown in parenthesis): HDF5 OMI/Aura Ozone (O3) Total Column 1-Orbit L2 Swath 13x24km (Level-2 OMTO3). I guess that you have other sort of data and GDAL can't georeference it and also gdalwarp will fail. The not-georeferenced tiff is probably the best that you can get but I may be wrong. — user30184, Apr 17 '19 at 19:59
If you look at the long list of ground control points from gdalinfo HDF5:"input.nc"://geophysical_data/Rrs_655 you'll see lots of these GCP[2676]: Id=, Info= (1560.5,7031.5) -> -88.2362823486328,43.7699966430664,0)? but also these GCP[2851]: Id=, Info= (260.5,7505.5) -> (-32767,-32767,0). They seem to present partly WGS84 coordinates and partly something totally different. — user30184, Apr 17 '19 at 20:43
The data file has latitude and longitude bands, so I suggest to follow https://gis.stackexchange.com/questions/154339/unable-to-warp-hdf5-files — AndreJ, Apr 19 '19 at 15:45

AndreJ · Accepted Answer · 2019-04-22T09:29:08.667

1

Gdalwarp stumbles over the nodata values in the latitude and longitude bands of the netcdf file. The related bug issue is fixed in GDAL 2.4.2: https://github.com/OSGeo/gdal/issues/1451

As a workaround, I did this:

Extract the desired band to a vrt:

 gdal_translate -of VRT HDF5:"input.nc"://geophysical_data/Rrs_655 -a_nodata -32767 input.vrt

Open the file in a text editor and remove all GCP that have coordinates of

X="-3.276700000000E+04" Y="-3.276700000000E+04"

Use gdalwarp WITHOUT the -geoloc parameter:

 gdalwarp -t_srs EPSG:4326 input.vrt input.tif

to get this result:

edited Apr 22 '19 at 09:29

answered Apr 21 '19 at 07:23

AndreJ

76,698
5
86
162

I was wondering if the above method valid for batch re-projection of netCDF files? – npal Apr 21 '19 at 18:26
Basically, yes. The vrt file is XML, and it should be possible to use an XML parser to eliminate the blocking GCP lines. But I can not help you on coding that in Python or C. – AndreJ Apr 21 '19 at 18:45
Thank you for the solution. I wanted to ask how did you choose those blocking GCP lines? I am trying to automate the process and was wondering if you could tell me how to certain which GCP lines to exclude from the code. – npal Apr 22 '19 at 13:25
You have to remove all GCP that have the nodata value -32767 for X and Y instead of useful degrees. For automation, you could read every GCP line, check for a valid range of +/-180, discard the GCPs outside of that, and write the others to a new vrt file. – AndreJ Apr 22 '19 at 16:46
Thank you AndreJ it works like a charm. The GDAL library has a bug which does not read the GCP and nodata values properly for HDF5 files. I have reported the issue and it has been fixed in the version 2.4.2. – npal Apr 22 '19 at 17:17
I am getting the image re-projected perfectly however the pixel values are a bit off. I have added the scale value as 0.000002 and offset value as 0.05. However, the pixels values still remain the same. They are in the range of -32767 to 0 and ideally they should between 0 to 1. – npal Apr 23 '19 at 15:07
Looks like you lost the nodata setting somewhere. – AndreJ Apr 23 '19 at 17:32
The commands that I am using are:
gdal_translate -of VRT HDF5:"input.nc"://geophysical_data/Rrs_655 -a_scale 0.000002 -a_offset 0.05 -a_nodata -32767 input6.vrt
gdalwarp -t_srs EPSG:4326 input6.vrt input5.tif
The nodata, scale and offset values are the same as that of the input.nc file. – npal Apr 23 '19 at 18:33
a_scale will not change pixel values; they still are INT16. I suggest to use gdal_calc to scale the values after reprojecting. And maybe cut the data to 25000/-30000 as mentioned in the metadata. – AndreJ Apr 24 '19 at 16:29
The command that I am using is: gdal_calc.py -A input1.tif --type='Float64' --outfile=result.tif --calc="A*0.000002+0.05". However, it seems to create a faulty image which I am unable to open in my software. – npal Apr 24 '19 at 18:17
What statistics has the resulting tif? Maybe you really need to clip invalid pixels. – AndreJ Apr 25 '19 at 05:40
The command that I used was: gdal_calc.py -A input7.tif --type='Float64' --outfile=result.tif --calc="A*0.000002+0.05". The files that are generated via the code are: 1. vrt file after removing the nodata GCPs : https://drive.google.com/open?id=1jV7CbY_W1qJFBL1ZowU73XIfjf4nE7Tj 2. The tif file generated after re-projection: https://drive.google.com/open?id=1D6HvCZz34hr0cARKL9gOOJokW8IrMcBU and 3. The result.tif generated after the gdal_calc is: https://drive.google.com/open?id=1D6HvCZz34hr0cARKL9gOOJokW8IrMcBU – npal Apr 25 '19 at 18:00
The third link is the same as the second. Anyway, the last tif should have values between 0.43 and 0.56. I don't know how your software handles nodata values. – AndreJ Apr 26 '19 at 06:19

Using GDALwarp for reprojecting netCDF file?

1 Answers1

Linked