3

My netcdf file coonsists of radar rainfall data in aeqd projection (azimuthal equidistant projection). Each file consists of 5952 timestamps(7.5 minutes resolution) of rainfall values for 1kmx1km grids.

I want to reproject from aeqd to UTM33 using gdalwarp in commandline as follows.

gdalwarp -t_srs "+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" -of netCDF C:/Radar/precip_aeqd-1000.201408.nc C:/Radar/RadarPrecip_UTM33proj.nc

I get error message that

Input file C:/Radar/precip_aeqd-1000.201408.nc has no raster bands.

Why do I get this error message?

I want to get output as reprojected netcdf file with utm33 1kmx1km grids.

Improve this question
12
  • 1
    You need to specify the subdataset to warp. See gis.stackexchange.com/questions/123107/… for the syntax. gdalinfo tells you the name of the subdatasets.
    – AndreJ
    Jun 3, 2017 at 10:49
  • By using gdalinfo, I found the subdatasets as follows. Subdataset_1 - equivalent_reflectivity_factor [5545x480x480] Subdataset_2 - lat Subdataset_3 - lon I am not clear how to specify the subdataset to warp.
    – Engkugan
    Jun 3, 2017 at 11:51
  • 1
    The srcfile part of gdalwarp has to be something like NETCDF:"path/to/file/filename":subdatasetname.
    – AndreJ
    Jun 3, 2017 at 12:30
  • I modified as you advised and run. I got the following warnings. Warning 1: NetCDF driver detected file type=5, but libnetcdf detected type=3. Warning: The target file has a 'nc' extension, which is normally used by the GMT, netCDF drivers, but the requested output driver is GTiff. Is it really what you want? Anyhow, I got the output file. However, when I open the file using netcdf4/Python, I get error message that file is not netcdf file format.
    – Engkugan
    Jun 6, 2017 at 9:17
  • Can you add the new command line to your question?
    – AndreJ
    Jun 6, 2017 at 9:51

1 Answer 1

Reset to default
1

I got it working in two steps for the first band:

gdal_translate -a_ullr -240500 240500 239500 -239500 -a_srs "+proj=aeqd +lat_0=59.458925 +lon_0=10.564472 +R=6371000 +x_0=0 +y_0=0 +units=m +no_defs" -a_nodata 9.96921e+036 -b 1 HDF5:"precip_aeqd-1000.201408.nc"://equivalent_reflectivity_factor precip-aeqd-1.tif
gdalwarp -overwrite -t_srs EPSG:25833 -tr 1000 1000 precip-aeqd.tif precip-utm-1.tif

This does not use the geoloc information in subdataset 2 and 3, but the corner coordinates from those subdatasets, converted to aeqd with

cs2cs +init=epsg:4326 +to +proj=aeqd +lat_0=59.458925 +lon_0=10.564472 +R=6371000 +x_0=0 +y_0=0 +units=m +no_defs corners.txt >aeqd.txt

and added 500m offset from mid cell to the edge. I omitted datum=WGS84 from the projection because this conflicts with the given radius.

This is my result, with the corner points from the subdatasets 2 and 3:

enter image description here

For all 5545 bands, you better use a vrt file:

gdal_translate -of VRT -a_ullr -240500 240500 239500 -239500 -a_srs "+proj=aeqd +lat_0=59.458925 +lon_0=10.564472 +R=6371000 +x_0=0 +y_0=0 +units=m +no_defs" -a_nodata 9.96921e+036 HDF5:"precip_aeqd-1000.201408.nc"://equivalent_reflectivity_factor precip-aeqd.vrt
gdalwarp -of GTiff -co COMPRESS=LZW -t_srs EPSG:25833 -tr 1000 1000 precip-aeqd.vrt precip-utm.tif

but the resulting tif is 2.4GB even with compression.

Improve this answer
8
  • Thank you very much. It is working well but I did some modifications.. When I use your aeqd projection text, output does not placed right location. I used aeqd projection from metadata and used it "+proj=aeqd +lat_0=59.458925 +lon_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs" and instead of 240000, I used 240500. It reprojected to right position.
    – Engkugan
    Jun 9, 2017 at 7:47
  • By using gdal_translate - of VRT, I got precip-aeqd.vrt and then used gdalwarp to the fianl output file precip-utm.tif. I get warning as follows. Warning 1: Metadata exceeding 32000 bytes cannot be written into GeoTIFF. Transferred to PAM instead.
    – Engkugan
    Jun 9, 2017 at 7:56
  • This warning is a restriction of the GeoTIFF format. Are you sure about +lon_0=0?
    – AndreJ
    Jun 9, 2017 at 7:57
  • As I asked initially, can I get the final output as netcdf file after reprojection?
    – Engkugan
    Jun 9, 2017 at 8:00
  • Are you sure about +lon_0=0? Yes,I think. I got the projection text from the file metadata. Further, I plotted the output in Arcmap. It placed on the right location.
    – Engkugan
    Jun 9, 2017 at 8:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.