Personally I don’t think this is a big problem, as long as there’s tiff that does the job of handling 8-bit rasters/images. Most of the netcdf grids I’ve been working with were elevation or bathymetry grids of floating point numbers, 8-bit unsigned integers weren’t of much use at all for these data. I am even wondering what is the use of 8-bit raster of Venus elevation with z values in range of 1 to 244. Mars elevation grid mentioned in the topic is beyond the range of 8-bit uint (values -7197 to +20834). But I am not experienced and have not been working with planetary datafiles 100Gb+ in size.
Hi Mikhail,
Let me clarify: These are not topography/elevation files, they are Synthetic Aperture Radar (SAR) image files, showing the radar brightness of the mapped regions of Venus. As such, z values in the 0-255 range are indeed the natural format, and they are not going to issue such files in any format other than unsigned 8-bit integer. And they shouldn’t issue them in any other format than unsigned 8-bit integer, because it keeps the file size at a sane level.
(I only used the Mars topography file as an example of a geoTIFF from the same group that was done right!)
And I agree, there shouldn’t be a problem if we can keep them as tiffs. But we can’t use the tiffs if the coordinate systems are messed up, and we can’t use most of the GMT/netCDF work-arounds that would normally work because of the unsigned 8-bit integer issue (amongst others). So we need to work with NASA/USGS to fix their issues at the source.
Sincerely,
Pat
sorry for persistence, but if the coordinate system is not correct, it does not matter whether its a 8-bit tiff or a floating point netcdf grid file format. And when the coordinate system is fixed, there is no problem working with 8-bit tiff using e.g. gdal tools.
I think we are saying the same thing, Mikhail! The coordinate system needs to be fixed, then everything will be fine.
In the long-term, that is.
Unfortunately, working with NASA/USGS to fix the issues will take time, probably many months. I am working with a Summer Intern on a project NOW, that only lasts 7 more weeks, and the problem is coming up with short-term GMT/netCDF work-arounds for the image problems. That’s what I’m lamenting above. We are now taking a different tack to solve our problem, forgive me if I don’t describe it in detail now, perhaps when we see some results next week…
Best regards,
Pat
It appears that gdalwarp can set georeferenced extent of the output grid to -180 -90 180 90 at the same time as changing coordinate system from Cartesian to geographical:
gdalwarp Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif \
-te -180 -90 180 90 -t_srs "+proj=lonlat +R=6051000 +no_defs" \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif -overwrite
One interesting line in its output is Creating output file that is 18775P x 9388L
Does this not literally mean pixel columns of Y and gridline rows of X?
...
ERROR 1: PROJ: eqc: Invalid latitude
Creating output file that is 18775P x 9388L.
Processing Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif [1/1] : 0Using internal nodata values (e.g. 0) for image Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif.
Copying nodata values from source Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif to destination Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif.
...10...20...30...40...50...60...70...80...90...100 - done.
the resulting ...-geo.tif appears to bear a correctly set coordinate grid, see gdalinfo and gmt grdinfo below.
gdalinfo Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif:
Driver: GTiff/GeoTIFF
Files: Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif
Size is 18775, 9388
Coordinate System is:
GEOGCRS["unknown",
DATUM["unknown",
ELLIPSOID["unknown",6051000,0,
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]],
PRIMEM["Reference meridian",0,
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.019174434087883,-0.019173412867490)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-180.0000000, 90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)
Upper Right ( 180.0000000, 90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)
Center ( 0.0000000, 0.0000000) ( 0d 0' 0.01"E, 0d 0' 0.01"N)
Band 1 Block=18775x1 Type=Byte, ColorInterp=Gray
NoData Value=0
gmt grdinfo Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif:
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif: Title: Grid imported via GDAL
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif: Command:
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif: Remark:
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif: Pixel node registration used [Geographic grid]
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif: Grid file format: gd = Import/export through GDAL
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif: x_min: -180 x_max: 180 x_inc: 0.0191744340879 name: x n_columns: 18775
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif: y_min: -90 y_max: 90 y_inc: 0.0191734128675 name: y n_rows: 9388
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif: v_min: 1.79769313486e+308 v_max: -1.79769313486e+308 name: z
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif: scale_factor: 1 add_offset: 0
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif: Default CPT:
+proj=longlat +R=6051000 +no_defs
Morover, the ...-geo.tif grid appears to be fully identical to the original tif (thanks to the default interpolator of gdalwarp I guess, which is nearest neighbor), as the difference grid ...-SUB.tif calculated using gmt grdmath -N is populated with zeroes, see gmt grdinfo output below:
gmt grdmath -N Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif \
SUB = Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-SUB.tif=gd:Gtiff
grdmath [WARNING]: x_inc does not divide 180; geographic boundary condition changed to natural.
gmt grdinfo Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-SUB.tif:
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-SUB.tif: Title: Grid imported via GDAL
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-SUB.tif: Command:
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-SUB.tif: Remark:
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-SUB.tif: Pixel node registration used [Geographic grid]
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-SUB.tif: Grid file format: gd = Import/export through GDAL
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-SUB.tif: x_min: -180 x_max: 180 x_inc: 0.0191744340879 name: x n_columns: 18775
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-SUB.tif: y_min: -90 y_max: 90 y_inc: 0.0191734128675 name: y n_rows: 9388
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-SUB.tif: v_min: 0 v_max: 0 name: z
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-SUB.tif: scale_factor: 1 add_offset: 0
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-SUB.tif: Default CPT:
+proj=longlat +R=6051000 +no_defs
18775*0.0191744340879 = 360; 0.0191744340879*(6051000*2*pi/360) = 2025.009549602
9388*0.0191734128675 = 180; 0.0191734128675*(6051000*2*pi/360) = 2024.901698645
plotting ...-geo.tif using grdimage still shows the same warning as gmt grdmath above:
gmt grdimage Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.tif -BWStr -Bxa -Bya --MAP_FRAME_TYPE=plain -png Venus-geo --MAP_FRAME_PEN=blue
grdimage [WARNING]: x_inc does not divide 180; geographic boundary condition changed to natural.
1- Still on the netCDF and char type. To be honest I don’t know the present state of it. I remember the many troubles I had along the years and that a char type didn’t even existed … till it was introduced several years ago. From Esteban’s post we still have issues with them and even gdalwarp when asked to do the inverse projection that Mikhail posted above but selecting .nc as output creates an image with type byte, that is a signed byte that ranges [-127 128].
2- The same gdalwarp prints lots of ERROR 1: PROJ: eqc: Invalid latitude which I think comes from the insurmountable problem that is the file has one more row that it should (that is why the the southern latitude is 30 arc sec below -90. What can be done bout this? Drop it? Disguise it by shortening the y_inc?
3-
The numbers make think so, but the GDAL model always create pixel-registered products. The answer, I think, lies in the fact that now x_inc != y_inc, that is, the disguising that I mentioned above.
I don’t think so. Please see below, gdalinfo of the .nc file produced by gdalwarp reports type as Byte, range as 0 to 255, and _Unsigned=true. The same does gmt grdinfo, v_min: 0 v_max: 255 name: GDAL Band Number 1
gdalwarp Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif \
-te -180 -90 180 90 -t_srs "+proj=lonlat +R=6051000 +no_defs" \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.nc -overwrite
gmt grdinfo Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.nc:
...
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.nc: Pixel node registration used [Geographic grid]
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.nc: Grid file format: nb = GMT netCDF format (8-bit integer), CF-1.7
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.nc: x_min: -180 x_max: 180 x_inc: 0.0191744340879 name: longitude n_columns: 18775
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.nc: y_min: -90 y_max: 90 y_inc: 0.0191734128675 name: latitude n_rows: 9388
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.nc: v_min: 0 v_max: 255 name: GDAL Band Number 1
...
but when plotted it can be and appears to be interpreted as signed:
gmt grdimage Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-geo.nc -BWEtr -Ba --MAP_FRAME_TYPE=plain -png Venus-geo-nc
What if there is not one more row but one column missing? I am returning back to my question regarding the first and the last column of a global grid, should these match? Considering gridline registration, the first column is centered on -180 degrees, while the last is centered at +180. This is essentially the same so the values should match, right? This is just an idea, maybe it is stupid, I don’t know how this is implemented in reality, and I am curious.
Then, assuming this “one column missing” idea is true brings us to the grid dimensions of 18776x9388 matching 2:1 ratio for a global geographical grid.
grid cell size is a square, side of 2024.9 meters.
Using my netCDF viewer (HDF Explorer) it tells me that the data type of the gdalwarp .nc file is NC_BYTE, and this is a 1 byte signed var
The other issue. GMT recognizes global grid registration grids with first col at -180 and last at +180 as well as the case where the last col is at +180-x_inc. First case, and this is from memory, even allows the case where first and last column are different (when they shouldn’t be), case in which, I think, both columns are averaged. But this is the GMT flexibility trying to recover from 3rd party liberties
Thanks for the explanations.
As a practical exercise I downloaded gmt’s 1 degree global earth relief data with grid registration, grid size 361x181. I extracted column 0 and column 360 using gdal_translate, the elevation values are identical, just longitude is 0 for col0 and 360 for col360. I understand that column 360 does not really need to be there explicitly.
This is not the case for pixel registered data. Raster 360x180, the first column is col0 at 0.5 degrees, the last column is col 359 is at 359.5 degrees and the last column is no duplicate of the first.
I decided to dig into the Venus grid the same way as I did with gmt 1-degree elevation grid.
First column of Venus data:
gdal_translate -srcwin 0 0 1 9388 \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif \
first_column.xyz
Last column of Venus data:
gdal_translate -srcwin 18774 0 1 9388 \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif \
last_column.xyz
Column before last:
gdal_translate -srcwin 18773 0 1 9388 \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif \
before_last_column.xyz
The last column is just filled with zeroes. Column before last is not identical to the first column.
so practically the raster x dimension appears to be 18774. More below.
First row of Venus data:
gdal_translate -srcwin 0 0 18774 1 \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif \
first_row.xyz
Last row of Venus data:
gdal_translate -srcwin 0 9387 18774 1 \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif \
last_row.xyz
Row before last:
gdal_translate -srcwin 0 9386 18774 1 \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif \
before_last_row.xyz
The first row contains data. Last row and the one before last are populated with zeroes. According to preview plots, there is large NaN area around the South pole so many rows of cells with zero values are expected.
Based on raster width of 18774, pixel width is 6051000 * 2 * 3.1415926 / 18774 = 2025.12 meters.
Seems reasonable to assume, based on 2:1 ratio, that possible real size of 2025x2025 m Venus raster is 18774x9387. This raster can be extracted from Venus grid by skipping that last column of zeroes and the last row of zeroes:
gdal_translate -srcwin 0 0 18774 9387 \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-18774x9387.tif
gdalwarp can be called to interpolate the 18774x9387 grid to geographical coordinates over the global region. Raster size -ts 18774 9387 must be specified explicitly, otherwise gdalwarp inherits coordinate system of the original 18775x9388 raster:
gdalwarp -ts 18774 9387 -te -180 -90 180 90 \
-t_srs "+proj=lonlat +R=6051000 +no_defs" \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-18774x9387.tif \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-18774x9387-geo.tif -overwrite
Or a new coordinate system can be assigned in place to the 18774x9387 grid (don’t know if this is the correct approach, but very fast and the plot of the resulting raster looks very similar to the output of gdalwarp):
gdal_edit.py -a_ullr -180 90 180 -90 \
-a_srs "+proj=lonlat +R=6051000 +no_defs" \
Venus_Magellan_C3-MDIR_Global_Mosaic_2025m-18774x9387.tif
So yes the grid looks buggy in many ways, as @Joaquim suspected.
The 2:1 argument is very compelling. I think you are right.
Thank you Mikhail and Joaquim for the ongoing discussion. In particular, the detailed analysis by Mikhail will be very helpful to me when I contact USGS Astrogeology to try to get them to fix their defective grid files.
I will note that I saw that row of zeros when I loaded the C3-MIDR into MATLAB a few days ago, and I was wondering how I was going to dispose of that if I re-wrote the grid. Our current work-around is to abandon the global grid and use a feature of the USGS website to request downloads of specific areas. We are looking at volcanic edifices and their surroundings, so we are typically looking at regions of 10° by 10° or so, and we avoid the grid edge issues that way.