search

OSGeo / PROJ-data

Repository for proj datum grids (for use by PROJ 7 or later)
Other
72 stars 33 forks source link

Updating GSIGEO2011 to include all the original data #50

Closed hernando closed 3 years ago

hernando commented 3 years ago

The initial commit was cropping the grid to the extent of EPSG:6695.

New the grid provides better results when compared to reference data:

Lat     Lon      geoid h (m) reference error (mm)
34.290  135.630  39.860329   39.8601    0.23
36.103  140.087  40.181506   40.1817   -0.19
43.217  143.129  30.638422   30.6389   -0.48
38.675  139.886  40.128380   40.1281    0.28
36.344  137.654  42.895627   42.8956    0.03
33.179  130.063  32.303562   32.3036   -0.04
39.801  141.322  41.886284   41.8862    0.08

The RMS is 0.637 mm (the original RMS was incorrectly reported as 0.9471, the correct value is 4.071).

hernando commented 3 years ago

I'm trying to regenerate the index with GDAL's master and I still have the problem of having a big diff (mostly due to slashes). I tried with a fresh checkout because I'm getting this problem even after doing make clean as Even suggested. I'm using Ubuntu 18.04, could that have anything to do with the problem?

rouault commented 3 years ago

I'm trying to regenerate the index with GDAL's master and I still have the problem of having a big diff (mostly due to slashes). I tried with a fresh checkout because I'm getting this problem even after doing make clean as Even suggested. I'm using Ubuntu 18.04, could that have anything to do with the problem?

You need GDAL 3.2 for the removal of the backslash escaping. I'll regenerate the index for you.

travis/expected_main.lst needs to be updated

rouault commented 3 years ago

./travis/check_new_grids.sh reports

The following warnings were found (the file will probably be read correctly by PROJ, but corrective action may be required):
 - Geotransform with a positive pixel height, that is a south-up image, is supported, but a unusual formulation

Could you run gdalwarp on the output of your initial conversion script, before optimizing the grid, to have a more conventional georeferencing that is with the TIFF image being north-up ?

kbevers commented 3 years ago

Any chance to get this PR finished today? I am releasing PROJ 7.2.0RC2 today, it would be great to include this update!

hernando commented 3 years ago

./travis/check_new_grids.sh reports

The following warnings were found (the file will probably be read correctly by PROJ, but corrective action may be required):
 - Geotransform with a positive pixel height, that is a south-up image, is supported, but a unusual formulation

Could you run gdalwarp on the output of your initial conversion script, before optimizing the grid, to have a more conventional georeferencing that is with the TIFF image being north-up ?

I have no clue about how to use gdalwarp for that. Instead I've reserved the order in which the data is processed to generate the first intermediate file. Now the grid looks correct to me.

Sorry, I still get backslashes with GDAL 3.2 (maybe because I combined it with PROJ 7.1.1?), please regenerate the index.

hernando commented 3 years ago

Sorry, I still get backslashes with GDAL 3.2 (maybe because I combined it with PROJ 7.1.1?), please regenerate the index.

Since it was a simple renaming I did it myself by hand to see if that fixes the build in Travis.

hernando commented 3 years ago

@rouault, could you take another look for merging?

hernando commented 3 years ago

Wait I have to double check the errors in the commit message, after the inversion they have changed

hernando commented 3 years ago

So, I have increased the number of significant digits of the resolution because I have the gut feeling that the actual grid is not using exact decimal precision for the location of the points. By doing that and leaving the grid inverted I actually got closer to the reference data I have:

Lat     Lon      geoid h (m) reference error (mm)
34.290  135.630  39.860122   39.8601   0.022
36.103  140.087  40.181748   40.1817   0.048
43.217  143.129  30.638869   30.6389   -0.031
38.675  139.886  40.128091   40.1281   -0.009
36.344  137.654  42.895648   42.8956   0.048
33.179  130.063  32.303559   32.3036   -0.041
39.801  141.322  41.886202   41.8862   0.002
RMS = 0.088 mm

Now, if I invert the order of the latitudes in the input I get this:

Lat     Lon      geoid h (m) reference error (mm)
34.290  135.630  39.864461   39.8601   4.361
36.103  140.087  40.177237   40.1817   -4.463
43.217  143.129  30.633345   30.6389   -5.555
38.675  139.886  40.132738   40.1281   4.638
36.344  137.654  42.895254   42.8956   -0.346
33.179  130.063  32.303632   32.3036   0.032
39.801  141.322  41.887449   41.8862   1.249
RMS = 9.643 mm

Almost 1 cm of RMS is too much in my opinion, don't you think?

gdalinfo tells me that the first grid has this extent

Upper Left  ( 119.9875000,  19.9916665) (119d59'15.00"E, 19d59'30.00"N)
Lower Left  ( 119.9875000,  50.0083332) (119d59'15.00"E, 50d 0'30.00"N)
Upper Right ( 150.0125000,  19.9916665) (150d 0'45.00"E, 19d59'30.00"N)
Lower Right ( 150.0125000,  50.0083332) (150d 0'45.00"E, 50d 0'30.00"N)
Center      ( 135.0000000,  34.9999998) (135d 0' 0.00"E, 35d 0' 0.00"N)

which is not good because it's upside down.

Whereas the inverted one has this:

Upper Left  ( 119.9875000,  50.0089335) (119d59'15.00"E, 50d 0'32.16"N)                                                 
Lower Left  ( 119.9875000,  19.9922668) (119d59'15.00"E, 19d59'32.16"N)                                                 
Upper Right ( 150.0125000,  50.0089335) (150d 0'45.00"E, 50d 0'32.16"N)                                                 
Lower Right ( 150.0125000,  19.9922668) (150d 0'45.00"E, 19d59'32.16"N)                                                 
Center      ( 135.0000000,  35.0006002) (135d 0' 0.00"E, 35d 0' 2.16"N)

which looks more correct, however GDAL has decided to adjust the latitudes very differently.

Could you tell me how to the the inversion with gdalwrap to see what result I get?

rouault commented 3 years ago

Could you tell me how to the the inversion with gdalwrap to see what result I get?

gdalwarp generates square pixels by default. You need to pass the explicit resolution you want with the -tr switch. On your original dataset, that is

gdalwarp ~/proj/PROJ-data/jp_gsi/jp_gsi_gsigeo2011_mainland.tif out.tif -overwrite -tr 0.025 0.016666700000000

(look at the Pixel Size reported by gdalinfo)

hernando commented 3 years ago

Looks good to me with the latest change. @hernando Is it ready from your perspective ?

Yes, I got distracted after the last push and forgot to let you know. It's interesting to see that adding more significant digits to the resolution of the latitude (in theory 1/60) gave better accuracy.

rouault commented 3 years ago

(I've refreshed files.geojson + index.html to reflect the larger extent, and updated dates) and the file is now on the CDN Quickly tested:

$ echo 35 140 0 | PROJ_NETWORK=ON PROJ_LIB=data src/cs2cs EPSG:6667 EPSG:6695
35.00   140.00 -34.94