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
Answer
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
No comments:
Post a Comment