Ticket #773 (closed defect: fixed)

Opened 5 years ago

Last modified 5 years ago

r.in.gdal unclean import of southern UTM zone

Reported by: hamish Owned by: grass-dev@…
Priority: major Milestone: 6.4.0
Component: Projections/Datums Version: svn-develbranch6
Keywords: r.in.gdal, utm Cc:
Platform: Unspecified CPU: Unspecified

Description

Hi,

after bringing in a GeoTiff? georef'd to UTM 48S with r.in.gdal I get a errors when trying to open the map:

WARNING: Raster map <Multichrom.red@user1> is in different zone
         (-48) than current region (48)
WARNING: Unable to open raster map <Multichrom.red@user1>
ERROR: Unable to open raster map <Multichrom.red>

location was created manually using the text based g.setproj.

r.in.gdal gave the "projection parameters appear to match" message.

G65> g.region -p
projection: 1 (UTM)
zone:       48
...

G65> g.proj -p
-PROJ_INFO-------------------------------------------------
name       : UTM
datum      : wgs84
towgs84    : 0.000,0.000,0.000
proj       : utm
ellps      : wgs84
a          : 6378137.0000000000
es         : 0.0066943800
f          : 298.2572235630
south      : defined
zone       : 48

GRASS 6.5svn, GDAL 1.5.2-3~bpo40+1 on Debian/etch.

gdalinfo:

...
Coordinate System is:
PROJCS["WGS 84 / UTM zone 48S",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
...

manually editing the $MAPSET/cellhd/ files to make the zone: line positive "fixes" it.

should the lib fn's automatically translate a cellhd negative zone number into PROJ_INFO's zone:+south: keys? cellhd can not use "south:" for utm zone as that's already used for the raster's southern bound. Should those two be allowed to get out of sync or is it less bad to just ignore the sign of the zone number in the raster's cellhd? If lib/gis/rd_cellhd.c is modified to allow negative zone numbers note that the uninit'd value of -1 would have to be adjusted to -255 or something outside the legal range.

?, Hamish

Change History

follow-up: ↓ 2   Changed 5 years ago by pkelly

Hi Hamish, How long ago was the location created? g.setproj had a bugfix on 7 October 2003 to create a negative zone number for southern hemisphere UTM locations (< http://freegis.org/cgi-bin/viewcvs.cgi/grass6/general/g.setproj/main.c.diff?r1=1.15&r2=1.16>). Are you saying that that doesn't work (i.e. the negative zone doesn't go into the WINDOW file) because libgis doesn't handle negative zone numbers? If so that would make this an incredibly long-running bug...

in reply to: ↑ 1   Changed 5 years ago by hamish

Replying to pkelly:

How long ago was the location created? g.setproj had a bugfix on 7 October 2003 to create a negative zone number for southern hemisphere UTM locations

It was created a few days ago with "grass65 -text" and a new location name, then defined by following the prompts.

Are you saying that that doesn't work (i.e. the negative zone doesn't go into the WINDOW file)

negative zone number doesn't make it into WIND, DEFAULT_WIND, or PROJ_INFO.

because libgis doesn't handle negative zone numbers? If so that would make this an incredibly long-running bug...

I've used southern hemisphere UTMs for ages without problems. I guess the "new" thing here is r.in.gdal imported data actually gets it right?

the PROJ_INFO file has

zone: (some positive number)
south: defined

so to me putting a negative zone in PROJ_INFO is just confusing+redundant, and so that part seems ok. (PROJ.4 uses +south to define a southern UTM btw)

The WIND and cellhd files can not have a "south: defined" tag as south: is used for the bound, so it should get a negative zone number there. r.in.gdal is doing that for cellhd files (but eg r.mapcalc does not for new maps), and WIND files are not getting it.

Part of the solution I think is to have G__.window.zone report a negative int for southern UTMs. That would mean that anything which writes a new raster or WIND file should do that too. (should opencell.c be using G_zone() instead?)

To avoid massive breakage from that we'd also have to improve the logic of G__open_cell_old() to test the abs() of the zone number and the south-ness of the projection (but how?).

Everything seems to run through the WIND file and not the PROJ_INFO file, so setting it in DEFAULT_WIND correctly seems to be the first step at fixing this properly.

which brings us back to g.setproj/main.c which is currently asking the 'South Hemisphere? [n]' question BEFORE the 'Enter Zone:' prompt, so when you enter the zone number it clobbers the cellhd.zone variable. The South question has to come after, somehow.

Hamish

  Changed 5 years ago by hamish

the following patch fixes the order of questions, and sets the WIND files correctly (and so new maps' cellhd too):

Index: general/g.setproj/proj-parms.table
===================================================================
--- general/g.setproj/proj-parms.table  (revision 39463)
+++ general/g.setproj/proj-parms.table  (working copy)
@@ -105,7 +105,7 @@
 UPS:Universal Polar Stereographic:SOUTH=ask,nodfl
 URM5:Urmaev V:ALPHA=ask,0.0;LAT0=ask,0.0;LON0=ask,20.0;NFACT=ask,1.0;QFACT=ask,1.0
 URMFPS:Urmaev Flat-Polar Sinusoidal:LAT0=ask,0.0;LON0=ask,20.0;NFACT=ask,1.0
-UTM:Universal Transverse Mercator (UTM):SOUTH=ask,nodfl;ZONE=ask,nodfl
+UTM:Universal Transverse Mercator (UTM):ZONE=ask,nodfl;SOUTH=ask,nodfl
 VANDG2:van der Grinten II:LAT0=ask,0.0;LON0=ask,20.0
 VANDG3:van der Grinten III:LAT0=ask,0.0;LON0=ask,20.0
 VANDG4:van der Grinten IV:LAT0=ask,0.0;LON0=ask,20.0

existing maps should be in existing mapsets/loc'ns so no bothered by this change.

ok to commit for trunk & 6.5?

Hamish

  Changed 5 years ago by hamish

  • status changed from new to closed
  • resolution set to fixed

committed proj-parms.table re-order in all branches; it doesn't seem to have broken anything. closing ticket.

Hamish

Note: See TracTickets for help on using tickets.