[GRASS-user] UTM 33x

mankoff · July 16, 2016, 11:42pm

Hi Grass List,

I have a newbie question about reprojecting data. I know how to import data if it has a coordinate system with an EPSG code, or even sometimes without. Now I have some data which is UTM 33X. I can import it just fine, with

grass70 -c -e file.grd ./grass/foo

But the region is then (from =g.region -p=),

projection: 0 (x,y)
zone: 0
north: 8562350
south: 8545750
west: 509450
east: 519550
nsres: 100
ewres: 100
rows: 166
cols: 101
cells: 16766

And if I'm in another projection, and I try to import with =r.proj=, it does not work:

WARNING: <PROJ_INFO> file not found for location <foo>
ERROR: Unable to get projection info of input map

The data provided told me this is UTM Zone 33X. There doesn't appear to be an EPSG code for that. Any advice will be much appreciated.

Thanks,

-k.

Micha_Silver · July 17, 2016, 5:45am

Hi

First, I have never heard of UTM 33X. Normally UTM zones are either north or south. So you might have UTM 33N, or UTM 33S. All UTM north zones that are based on the WGS84 datum are numbered: 326zz, where zz is the zone. So if your data is in UTM 33N, the EPSG code will be 32633. All UTM south projections (again based on WGS84) are numbered 327zz, with the zz=zone again. Next, g.region does not affect the projection. Only the region (extents and resolution) within the current LOCATION’s projection. The r.proj module does not import. Rather it does reprojection between one CRS and another. So: I think you should setup a GRASS LOCATION and MAPSET (using the wizard) based on the EPSG code 32633 (or 32733 if your data is UTM 33S), then, within a GRASS session in that LOCATION just do:

···

On 17/07/2016 02:42, Ken Mankoff wrote:

Hi Grass List,

I have a newbie question about reprojecting data. I know how to import data if it has a coordinate system with an EPSG code, or even sometimes without. Now I have some data which is UTM 33X. I can import it just fine, with

grass70 -c -e file.grd ./grass/foo

But the region is then (from =g.region -p=),

projection: 0 (x,y)
zone:       0
north:      8562350
south:      8545750
west:       509450
east:       519550
nsres:      100
ewres:      100
rows:       166
cols:       101
cells:      16766

And if I'm in another projection, and I try to import with =r.proj=, it does not work:

WARNING: <PROJ_INFO> file not found for location <foo>
ERROR: Unable to get projection info of input map

The data provided told me this is UTM Zone 33X. There doesn't appear to be an EPSG code for that. Any advice will be much appreciated.

r.in.gdal -c input=file.grd output=file

Thanks,

  -k.
_______________________________________________
grass-user mailing list
[grass-user@lists.osgeo.org](mailto:grass-user@lists.osgeo.org)
[http://lists.osgeo.org/mailman/listinfo/grass-user](http://lists.osgeo.org/mailman/listinfo/grass-user)

-- 
Micha Silver
Arava Drainage Authority
+972-523-665918

mankoff · July 17, 2016, 11:06am

Hi Micha,

Sorry if my terminology was incorrect and confused the issue.

On 2016-07-17 at 05:45, Micha Silver <micha@arava.co.il> wrote:

First, I have never heard of UTM 33X. Normally UTM zones are either
north or south. So you might have UTM 33N, or UTM 33S. All UTM north
zones that are based on the WGS84 datum are numbered: 326zz, where zz
is the zone. So if your data is in UTM 33N, the EPSG code will be
32633. All UTM south projections (again based on WGS84) are numbered
327zz, with the zz=zone again.

33X exists. It covers Svalbard.

See https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Exceptions

> And if I'm in another projection, and I try to import with =r.proj=,
> it does not work:

The r.proj module does not import. Rather it does reprojection between
one CRS and another.

Yes, I mean reproject, although I think of it as importing from an existing location to a new location.

So: I think you should setup a GRASS LOCATION and MAPSET (using the
wizard) based on the EPSG code 32633 (or 32733 if your data is UTM
33S), then, within a GRASS session in that LOCATION just do: r.in.gdal
-c input=file.grd output=file

I've let grass define the location using the "-c" flag to the grass command, then imported it with r.in.gdal. But because it is 33X, the location isn't fully defined, and I cannot reproject.

If I follow your advice, it does not work:

$ grass70 -c epsg:32633 tmp
GRASS 7.0.3 (tmp): > r.in.gdal input=file.grid output=file
ERROR: Projection of dataset does not appear to match current location.

-k.

hellik · July 17, 2016, 11:36am

33X exists. It covers Svalbard.

See

https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Exceptions

Grids & Datums: THE KINGDOM OF NORWAY

http://www.asprs.org/a/resources/grids/10-99-norway.pdf

maybe it helps to define a GRASS location for this area for your own.

also have a look to www.epsg.io to find maybe an ESPSG code for this UTM
zone.

-----
best regards
Helmut
--
View this message in context: http://osgeo-org.1560.x6.nabble.com/UTM-33x-tp5276714p5276744.html
Sent from the Grass - Users mailing list archive at Nabble.com.

mankoff · July 17, 2016, 12:10pm

Hi Helmut,

On 2016-07-17 at 11:36, Helmut Kudrnovsky <hellik@web.de> wrote:

33X exists. It covers Svalbard.

Grids & Datums: THE KINGDOM OF NORWAY

http://www.asprs.org/a/resources/grids/10-99-norway.pdf

After a quick skim of that article, I'm not sure why you sent it. It doesn't mention 33X, and only,

The coordinate systems of Svalbard include the Thumb Point Datum of
1948 origin where: Theta_0 79° 03’ 58.97” North, A_0 2° 48’ 15.36”
East of Greenwich, and is referenced to the Bessel 1841 ellipsoid.

But the data I have is referenced to WGS 84. It is UTM 33X, which is 72°N to 84°N, and is 12° wide, from 9° to 21°.

also have a look to www.epsg.io to find maybe an ESPSG code for this UTM
zone.

There is none that I can find.

maybe it helps to define a GRASS location for this area for your own.

I'm not sure what you mean by "define a GRASS location for this area for your own". I think I've done that (or more correctly GRASS has done that):

grass70 -text -c -e file.grid ./grass/foo
grass70 -text ./grass/foo/PERMANENT
r.import input=file.grd output=f

I have a grass location. The data is imported. I can do analysis on it and visualize it. But being a newbie to grass, I don't know how to get it out of this location and into another location for analysis that requires other data.

Because 33X is well defined (see my description above, based on this map http://www.dmap.co.uk/utmworld.htm ), I assume it is possible to export or reproject it to a different location, but being a newbie I don't know how to do this in GRASS, hence these emails.

-k.

hellik · July 17, 2016, 12:18pm

http://www.asprs.org/a/resources/grids/10-99-norway.pdf

After a quick skim of that article, I'm not sure why you sent it. It

doesn't mention 33X, and only,

in these documents there are sometimes nice informations about special cases
of projections, grids etx in some regions of the world.

maybe it helps to define a GRASS location for this area for your own.

I'm not sure what you mean by "define a GRASS location for this area for

your own".

you can define a location by define all needed parameters manually, have a
look at the GRASS startup wizzard

I think I've done >that (or more correctly GRASS has done that):

grass70 -text -c -e file.grid ./grass/foo
grass70 -text ./grass/foo/PERMANENT
r.import input=file.grd output=f

please post the output of g.proj -p
(https://grass.osgeo.org/grass73/manuals/g.proj.html)

and also post the output of: gdalinfo file.grid

or try: testepsg file.grid (see attached the file
gdal_testepsg_espsg_32633.txt
<http://osgeo-org.1560.x6.nabble.com/file/n5276749/gdal_testepsg_espsg_32633.txt>
of the example testepsg EPSG:32633 to get an idea which parameters are
needed to define a projection).

But being a newbie to grass, I don't know how to get it out of this

location and into another location for >analysis that requires other data.

the outputs of g.proj -p and/or gdalinfo file.grid are needed to see which
parameters may be missing for a successfull reprojection into another
location.

-----
best regards
Helmut
--
View this message in context: http://osgeo-org.1560.x6.nabble.com/UTM-33x-tp5276714p5276749.html
Sent from the Grass - Users mailing list archive at Nabble.com.

mankoff · July 17, 2016, 12:51pm

On 2016-07-17 at 12:18, Helmut Kudrnovsky <hellik@web.de> wrote:

you can define a location by define all needed parameters manually,
have a look at the GRASS startup wizard

1) grass70 -gui
2) Select GRASS location
3) [new]
4) [Next>]
5) (o) Create a new Cartesian coordinate system (XY)
6) [Next>]
7) [Finish]
8) Do you want to set the default region extents now? [Yes]
9) Set 72,84,9,21 as the S,N,W,E boundaries
10) [Set Region]

GRASS freezes. When I kill the process, terminal prints:

ERROR: Invalid return code from GUI startup script.
Please advise GRASS developers of this error.
Exiting...

please post the output of g.proj -p
(https://grass.osgeo.org/grass73/manuals/g.proj.html)

XY location (unprojected)

(that is from =g.proj -p= run in the GRASS location where the GRASS CLI creates the location via the "-c" flag and the file I then import.

and also post the output of: gdalinfo file.grid

Driver: GS7BG/Golden Software 7 Binary Grid (.grd)
Files: file.grd
Size is 101, 166
Coordinate System is `'
Origin = (509450.000000000000000,8562350.000000000000000)
Pixel Size = (100.000000000000000,-100.000000000000000)
Corner Coordinates:
Upper Left ( 509450.000, 8562350.000)
Lower Left ( 509450.000, 8545750.000)
Upper Right ( 519550.000, 8562350.000)
Lower Right ( 519550.000, 8545750.000)
Center ( 514500.000, 8554050.000)
Band 1 Block=101x1 Type=Float64, ColorInterp=Undefined
  Min=-2.316 Max=385.699
  NoData Value=1.70141000000000007e+38

or try: testepsg file.grid (see attached the file
gdal_testepsg_espsg_32633.txt

ERROR 1: Error occured translating file.grd.

-k.

hellik · July 17, 2016, 12:46pm

please post the output of g.proj -p
(https://grass.osgeo.org/grass73/manuals/g.proj.html)

XY location (unprojected)

(that is from =g.proj -p= run in the GRASS location where the GRASS CLI

creates the location via the "-c" flag and the file I then import.

and also post the output of: gdalinfo file.grid

Driver: GS7BG/Golden Software 7 Binary Grid (.grd)
Files: file.grd
Size is 101, 166
Coordinate System is `'
Origin = (509450.000000000000000,8562350.000000000000000)
Pixel Size = (100.000000000000000,-100.000000000000000)
Corner Coordinates:
Upper Left ( 509450.000, 8562350.000)
Lower Left ( 509450.000, 8545750.000)
Upper Right ( 519550.000, 8562350.000)
Lower Right ( 519550.000, 8545750.000)
Center ( 514500.000, 8554050.000)
Band 1 Block=101x1 Type=Float64, ColorInterp=Undefined
Min=-2.316 Max=385.699
NoData Value=1.70141000000000007e+38

or try: testepsg file.grid (see attached the file
gdal_testepsg_espsg_32633.txt

ERROR 1: Error occured translating file.grd.

and from your original report in this thread:

WARNING: <PROJ_INFO> file not found for location <foo>
ERROR: Unable to get projection info of input map

as you can see your data has no projection information, thus imported to an
xy location (that means a location without any projection info).

GRASS GIS (or any other GIS) can't re-project such data to another
projection.

some work arounds:

- ask the data provider for data _with_ projection informations
- have a look at e.g. [1], assign projection information to your raster by
e.g. gdal_translate

HTH

[1] https://lists.osgeo.org/pipermail/gdal-dev/2009-November/022607.html

-----
best regards
Helmut
--
View this message in context: http://osgeo-org.1560.x6.nabble.com/UTM-33x-tp5276714p5276752.html
Sent from the Grass - Users mailing list archive at Nabble.com.

Micha_Silver · July 17, 2016, 1:17pm

Hi Ken

Well, shame on me for missing that. I like to take part in these mail lists because I learn something new all the time! Helmut already pointed to trying gdalinfo on the file.grid. If that does not return a full projection description then either get the full projection info from the data provider, or try to manually setup a GRASS location with the parameters for UTM 33X, then import the original into that correctly defined location with the -o option to r.in.gdal . At that point you will be able to re-project to other projections.

···

On 07/17/2016 02:06 PM, Ken Mankoff wrote:

Hi Micha,

Sorry if my terminology was incorrect and confused the issue.

On 2016-07-17 at 05:45, Micha Silver [<micha@arava.co.il>](mailto:micha@arava.co.il) wrote:

First, I have never heard of UTM 33X. Normally UTM zones are either
north or south. So you might have UTM 33N, or UTM 33S. All UTM north
zones that are based on the WGS84 datum are numbered: 326zz, where zz
is the zone. So if your data is in UTM 33N, the EPSG code will be
32633. All UTM south projections (again based on WGS84) are numbered
327zz, with the zz=zone again.

33X exists. It covers Svalbard.

See [https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Exceptions](https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Exceptions)

And if I'm in another projection, and I try to import with =r.proj=,
it does not work:

The r.proj module does not import. Rather it does reprojection between
one CRS and another.

Yes, I mean reproject, although I think of it as importing from an existing location to a new location.

So: I think you should setup a GRASS LOCATION and MAPSET (using the
wizard) based on the EPSG code 32633 (or 32733 if your data is UTM
33S), then, within a GRASS session in that LOCATION just do: r.in.gdal
-c input=file.grd output=file

I've let grass define the location using the "-c" flag to the grass command, then imported it with r.in.gdal. But because it is 33X, the location isn't fully defined, and I cannot reproject.

If I follow your advice, it does not work:

$ grass70 -c epsg:32633 tmp
GRASS 7.0.3 (tmp): > r.in.gdal input=file.grid  output=file
ERROR: Projection of dataset does not appear to match current location.

  -k.

mankoff · July 17, 2016, 1:32pm

Hi Helmut,

Thanks for the advice so far.

On 2016-07-17 at 12:46, Helmut Kudrnovsky <hellik@web.de> wrote:

as you can see your data has no projection information, thus imported
to an xy location (that means a location without any projection info).

GRASS GIS (or any other GIS) can't re-project such data to another
projection.

If I could get GRASS to define the lat,lon boundaries as I tried to do in my last email, but without crashing, wouldn't that solve this?

And perhaps I'm still confused, but the projection is *known*, even if it isn't encoded in this file. If this were known as 33N instead of 33X, I could define a location with an EPSG code, and then import with "r.in.gdal -o", where "-o" is "Override projection (use location's projection)". It seems like I should be able to do something similar, but I may need to define the projection with a WKT file or proj4 code which can support 33X.

-k.

hellik · July 17, 2016, 1:36pm

If I could get GRASS to define the lat,lon boundaries as I tried to do in my

last email, but without >crashing, wouldn't that solve this?

mybe yes

And perhaps I'm still confused, but the projection is *known*, even if it

isn't encoded in this file.

that's the point: the file itself is missing the projection information, see
your gdalinfo output:

"Coordinate System is `'"

and this information is needed for accurate re-projection for GIS.

If this were known as 33N instead of 33X, I could define a location with an

EPSG code, and then import >with "r.in.gdal -o", where "-o" is "Override
projection (use location's projection)".

that would be the first what I may do (as mentioned in the GDAL ML thread):
try the data to import into UTM 33N location and see if there are
distortions or any other problems at a first step.

another idea:

ask in the proj-ML (http://lists.maptools.org/mailman/listinfo/proj) and
maybe there metacrs (https://lists.osgeo.org/mailman/listinfo/metacrs).

the former is known for good advices for defining projections.

HTH

-----
best regards
Helmut
--
View this message in context: http://osgeo-org.1560.x6.nabble.com/UTM-33x-tp5276714p5276758.html
Sent from the Grass - Users mailing list archive at Nabble.com.

mankoff · July 17, 2016, 2:48pm

Hi Micha and Helmut,

Thanks again for your help so far, and sorry to everyone else for filling your mailboxes on a Sunday morning.

It seems that Micha suggested the fix several emails ago (as did Helmut), which is to try to import into 33N and see what happens. I replied:

If I follow your advice, it does not work:
$ grass70 -c epsg:32633 tmp
GRASS 7.0.3 (tmp): > r.in.gdal input=file.grid output=file
ERROR: Projection of dataset does not appear to match current location.

But if you tell GRASS to ignore the data projection, with "r.in.gdal -o", then it does appear to work (especially after dealing with other issues, such as re-setting the region with "g.region raster=file" after importing).

Thank you, Micha and Helmut, for helping me figure this out,

-k.

Blumentrath_Stefan1 · July 18, 2016, 7:34am

Hei,

Even if the issue is solved, some additional info might shed some more light on this.

The X in 33X is specific for the Military Grid Reference System (http://mgrs-data.org/), which is still frequently used in e.g. some botanical journals.

There is no EPSG code for MGRS since there are several exceptions for grid cell sizes where cells like 33X or 32V were expanded, while other, neighboring cells were shrunk.

For practical purposes, 33X can be treated as WGS84/UTM 33N represented by EPSG:32633 (as e.g. 32V equals to WGS84/UTM 32N).

Kind regards,
Stefan

-----Original Message-----
From: grass-user [mailto:grass-user-bounces@lists.osgeo.org] On Behalf Of Ken Mankoff
Sent: 17. juli 2016 16:48
To: Micha Silver <micha@arava.co.il>
Cc: GRASS users list <grass-user@lists.osgeo.org>
Subject: Re: [GRASS-user] UTM 33x