[GRASSLIST:10027] projection problems

GRASS GRASS user
1

It seems that I could use some very basic explanation concerning locations and projections. I am trying to read several shape files into a single location. However, I am getting a projection error on some of them.

g.region -p for the location is:

GRASS 6.1.cvs (minnesota_utm):~ > g.region -p
projection: 1 (UTM)
zone: 15
datum: nad83
ellipsoid: grs80
north: 5472415.17999963
south: 4816310.42
west: 189788.954
east: 761663.27300056
nsres: 0.99999963
ewres: 1.00000056
rows: 656105
cols: 571874

and v.info for the state outline shape file (which read in fine with v.in.ogr) is:

GRASS 6.1.cvs (minnesota_utm):~ > v.info mn_state
+----------------------------------------------------------------------------+
| Layer: mn_state Organization: |
| Mapset: pls Source Date: |
| Location: minnesota_utm Name of creator: |
| Database: /Volumes/disk1/home1/kwythers/grassdata |
| Title: |
| Map Scale: 1:1 |
| Map format: native |
|----------------------------------------------------------------------------|
| Type of Map: Vector (level: 2) |
| |
| Number of points: 0 Number of areas: 1 |
| Number of lines: 0 Number of islands: 1 |
| Number of boundaries: 1 Number of faces: 0 |
| Number of centroids: 1 Number of kernels: 0 |
| |
| Map is 3D: 0 |
| Number of dblinks: 1 |
| |
| Projection: UTM (zone 0) |
| N: 5472414.182 S: 4816310.425 |
| E: 761662.273 W: 189788.954 |
| B: 0.000 T: 0.000 |
| |
| Digitize threshold: 0.00000 |
| Comments: |
| |
+----------------------------------------------------------------------------+

Although I do notice that the projection reads UTM (zone 0) which seems odd to me... it should be zone 15.

Now for the really fun part...

If I try and read another shape file into the location, I get the projection error about not matching the current location:

GRASS 6.1.cvs (minnesota_utm):~ > v.in.ogr dsn=~/projects/mn_dnr/shape_files/pls_corners_btrees/county/cook/ output=cook_test
ERROR: Projection of dataset does not appear to match current location.

        LOCATION PROJ_INFO is:
        name: UTM
        datum: nad83
        towgs84: 0.000,0.000,0.000
        proj: utm
        ellps: grs80
        a: 6378137.0000000000
        es: 0.0066943800
        f: 298.2572221010
        zone: 15

        cellhd.proj = 0 (unreferenced/unknown)

        You can use the -o flag to v.in.ogr to override this check.
        Consider to generate a new location with 'location' parameter from
        input data set.

If I use the -o flag (or location=some_new_location) the file comes in as projection (0) x,y

I asked one of our ARC/GIS folks downstairs to see if he could import the shape files and lay them on top of each other and he did not have any trouble. I am wondering if I have set something very basic up wrong.

Any suggestions for diagnosing this problem and getting these shape files to match up would be much appreciated.

Thanks

2

On Fri, Jan 27, 2006 at 11:52:04AM -0600, we recorded a bogon-computron collision of the <kwythers@umn.edu> flavor, containing:

It seems that I could use some very basic explanation concerning
locations and projections. I am trying to read several shape files
into a single location. However, I am getting a projection error on
some of them.

[...]

If I try and read another shape file into the location, I get the
projection error about not matching the current location:

GRASS 6.1.cvs (minnesota_utm):~ > v.in.ogr dsn=~/projects/mn_dnr/
shape_files/pls_corners_btrees/county/cook/ output=cook_test
ERROR: Projection of dataset does not appear to match current location.

       LOCATION PROJ_INFO is:
       name: UTM
       datum: nad83
       towgs84: 0.000,0.000,0.000
       proj: utm
       ellps: grs80
       a: 6378137.0000000000
       es: 0.0066943800
       f: 298.2572221010
       zone: 15

       cellhd.proj = 0 (unreferenced/unknown)

       You can use the -o flag to v.in.ogr to override this check.
       Consider to generate a new location with 'location' parameter
from
       input data set.

I asked one of our ARC/GIS folks downstairs to see if he could import
the shape files and lay them on top of each other and he did not have
any trouble. I am wondering if I have set something very basic up wrong.

Do your shapefiles have ".prj" files to go with them? My guess is they do
not.

If not, there's no way for v.in.ogr to know what projection they're in, and it
will always say that there's a projection mismatch. .prj files hold "WKT"
specifications of the coordinate system/projection of the associated shapefile. They are not always distributed along with shapefiles --- a fact that probably
casts a dark shadow on the souls of members of GIS departments that distribute
those files and leave folks to guess at what type of coordinates lie inside.

Odds are good that your ArcGIS friends are able to import them because ArcGIS
is simply reading in the shape file and *assuming* it's in the right projection
for the current view. That's what "v.in.ogr -o" is supposed to do, too.
As long as that assumption's good, the data will lay on top of each other
properly.

If I use the -o flag (or location=some_new_location) the file comes
in as projection (0) x,y

If you're sure the shapefile is in UTM zone 15 coordinates with NAD83 datum,
you should be able to use v.in.ogr with -o but without "location=newlocation"
and grass should ignore the projection mismatch and import the shapes
with the assumption that the coordinates in the shapefile match the coordinate
system of your location. It should not be possible to get a "projection (0)
x,y" when importing a shapefile into an existing UTM location, but if you use
"location=newlocation" then that's just what it'll do because it has no
information about what projection was used to create the shapefile.

Any suggestions for diagnosing this problem and getting these shape
files to match up would be much appreciated.

If you *know for certain* that your shapefiles definitely contain shapes
defined in terms of NAD83 UTM zone 15 coordinates, you should be able to
import them into your NAD83 UTM zone 15 GRASS location with v.in.ogr -o.
Or, you can give v.in.ogr a hand by providing it with a .prj file to go
along with each of the shapefiles. The WKT for UTM zone 15 NAD83 in the
northern hemisphere is:

PROJCS["NAD83 / UTM zone 15N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-93],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1]]

If you put this single line of text in a file called "foo.prj" then it will
be used to define the coordinate system when you import foo.shp. From then
on GRASS would be able to tell immediately whether the shapefile matched the
projection of the current location, and other ogr tools (like ogr2ogr) could
be used pretty simply to reproject or datum-shift the files.

It ought to be illegal to distribute shapefile sets without .prj files, but
many ESRI houses do just that.

Naturally, it would be wrong to use this .prj file with a shapefile that was
defined in some other coordinates --- NAD27 UTM zone 15, for example. But
that's another story.

--
Tom Russo KM5VY SAR502 DM64ux http://www.swcp.com/~russo/
Tijeras, NM QRPL#1592 K2#398 SOC#236 AHTB#1
"The only thing you can do easily is be wrong, and that's hardly
  worth the effort." -- Norton Juster

3

Thank you very much for the reply Tom. My responses are in line.

On Jan 27, 2006, at 1:58 PM, Tom Russo wrote:

Do your shapefiles have ".prj" files to go with them? My guess is they do
not.

You are correct.

If not, there's no way for v.in.ogr to know what projection they're in, and it
will always say that there's a projection mismatch. .prj files hold "WKT"
specifications of the coordinate system/projection of the associated shapefile. They are not always distributed along with shapefiles --- a fact that probably
casts a dark shadow on the souls of members of GIS departments that distribute
those files and leave folks to guess at what type of coordinates lie inside.

Odds are good that your ArcGIS friends are able to import them because ArcGIS
is simply reading in the shape file and *assuming* it's in the right projection
for the current view. That's what "v.in.ogr -o" is supposed to do, too.
As long as that assumption's good, the data will lay on top of each other
properly.

If I use the -o flag (or location=some_new_location) the file comes
in as projection (0) x,y

If you're sure the shapefile is in UTM zone 15 coordinates with NAD83 datum,
you should be able to use v.in.ogr with -o but without "location=newlocation"
and grass should ignore the projection mismatch and import the shapes
with the assumption that the coordinates in the shapefile match the coordinate
system of your location. It should not be possible to get a "projection (0)
x,y" when importing a shapefile into an existing UTM location, but if you use
"location=newlocation" then that's just what it'll do because it has no
information about what projection was used to create the shapefile.

Any suggestions for diagnosing this problem and getting these shape
files to match up would be much appreciated.

If you *know for certain* that your shapefiles definitely contain shapes
defined in terms of NAD83 UTM zone 15 coordinates, you should be able to
import them into your NAD83 UTM zone 15 GRASS location with v.in.ogr -o.

I went back and re-imported with v.in.ogr -o. That worked in that v.info output makes sense now for both the state outline and the vector points files (those were the troublesome ones). And d.vect can plot both the state outline and the points on top of each other, and the appearance... to my eye, makes sense. Which is a good thing.

However, I notice in the points file (I'm not sure what else to call it) which is a grass vector containing site information in a postgresql table, of which two columns 'x' and 'y' contain the coordinates for plotting each site, do not make sense for Minnesota UTM coordinates. For example, v.info on the file gives:

  Projection: UTM (zone 0)
| N: 5470474.000 S: 4816355.500
| E: 760498.938 W: 189916.844
| B: 0.000 T: 0.000

As you can see the North and South dimension show UTM coordinates are in the 4 to 5 million range. However checking the 'x' and 'y' columns in the postgresql table the values are an order on magnitude too low. (ie x=441988 and y=508562).

Is this due to the use of the -o flag? In that you were saying that the use of -o allows v.in.ogr to assume that the projection is ok and forces the points to plot within the projection's dimensions?

If that is the case, then that explains why I could not plot the points file I created from a sql cross join command. The results of which pulled each of the 4 species out of a single row (in a 250k row file), and gave each species its own row (creating a 1 million row file). Then I read the new table back into grass with v.in.db (which has no -o option... that I am aware of). That new vector then would not plot with the other data in the region.

Now I have not re-done the sql cross join command on the vector points file yet. It is running right now (takes a while). If I understand your explanation Tom, I do not expect it to work any better than it did last time.

I guess the big question right now is why would the shape file I am trying to read have metadata describing UTM zone 15 nad83 bla... bla... bla..., and yet be full on point data that could not possible be UTM coordinates for Minnesota? Then what can I do about it? Would using your .prj info below solve this?

Thanks again for taking the time to write out such a lengthy explanation.

Kirk

Or, you can give v.in.ogr a hand by providing it with a .prj file to go
along with each of the shapefiles. The WKT for UTM zone 15 NAD83 in the
northern hemisphere is:

PROJCS["NAD83 / UTM zone 15N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-93],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1]]

If you put this single line of text in a file called "foo.prj" then it will
be used to define the coordinate system when you import foo.shp. From then
on GRASS would be able to tell immediately whether the shapefile matched the
projection of the current location, and other ogr tools (like ogr2ogr) could
be used pretty simply to reproject or datum-shift the files.

It ought to be illegal to distribute shapefile sets without .prj files, but
many ESRI houses do just that.

Naturally, it would be wrong to use this .prj file with a shapefile that was
defined in some other coordinates --- NAD27 UTM zone 15, for example. But
that's another story.

--
Tom Russo KM5VY SAR502 DM64ux http://www.swcp.com/~russo/
Tijeras, NM QRPL#1592 K2#398 SOC#236 AHTB#1
"The only thing you can do easily is be wrong, and that's hardly
  worth the effort." -- Norton Juster

4

Tom, before you spend too much time on this...

I think DNR data came with a screwy y value in their table. I opened the original DNR table, and as I've said before, the values are there that are all less than 1 million - which can't be correct. I added x and y coords, named something like xadded, yadded (bad short term memory), and then calculated the difference, saved in ydiff. It is a number that is something like 4,700,206, plus or minus 0.5.

I have no idea why MN DNR did this. Do you see the same thing? But I am beginning to think it is a not grass or postgresql issue...

On Jan 30, 2006, at 10:05 AM, Tom Russo wrote:

5

On Jan 30, 2006, at 10:05 AM, Tom Russo wrote:

Your v.info output shows what looks like valid UTM zone 15 data for
Minnesota. I can't tell you why your database has strange values in it.
I have not made much use of grass with postgres, so I don't have much
experience to understand that process. Whenever I do a v.in.ogr it doesn't
make a postgres table for me.

I started using postgres to handle attributes so that I could have access to sql querries. Here is what I do: First make location and mapset, then use db.connect to connect to postgresql database (the database needs to exist in postgres), then use v.in.ogr to read in shape files. The tables all get made automatically in the postgres database (just like they would in dbf is that was the database connection with db.connect).

Now here is what I think is going on with these damn files... tell me if this makes sense to you. The coordinates in the .shp file from the state are correct. This is why if I read the files in with v.in.ogr, they can be displayed with d.vect on top of an outline of the state is all looks fine. However, the y coordinates in the attribute file (which come from the dbf file that came with the .shp file from the state), are all messed up by some 4,700,206, plus or minus 0.5 (for god knows what reason).

Here is what I'm wondering... if the .shp file has the correct geography coordinates, is there some command that out write out those coordinates? If so, I could substitute the good y coordinates for the bad y coordinates in the attribute table. I assume if I use v.out.ogr, that just writes the attributes from the attribute table (including the weird y coordinates). true?

Kirk

Could you possibly put your shapefiles (with metadata) on a web site or ftp
site so that I can download them and take a look?

--
Tom Russo KM5VY SAR502 DM64ux http://www.swcp.com/~russo/
Tijeras, NM QRPL#1592 K2#398 SOC#236 AHTB#1
"The only thing you can do easily is be wrong, and that's hardly
  worth the effort." -- Norton Juster

6

On Mon, Jan 30, 2006 at 11:40:06AM -0600, we recorded a bogon-computron collision of the <kwythers@umn.edu> flavor, containing:

Tom, before you spend too much time on this...

I think DNR data came with a screwy y value in their table. I opened
the original DNR table, and as I've said before, the values are there
that are all less than 1 million - which can't be correct. I added x
and y coords, named something like xadded, yadded (bad short term
memory), and then calculated the difference, saved in ydiff. It is a
number that is something like 4,700,206, plus or minus 0.5.

Yes, it does sound like your DNR table is not in the same coordinates as
your other data. Without more information, of course, I couldn't do more
than guess at what's going on. It could even be that your DNR table is in
Minnesota state plane coordinate system rather than UTM. I think that
in any of the three Minnesota SPCS coordinate systems the northing would
be an order of magnitude smaller than they'd be in UTM. The only way to
resolve that question would be to check, double check, and triple check with
the office that provided the data.

I have no idea why MN DNR did this. Do you see the same thing? But I
am beginning to think it is a not grass or postgresql issue...

It does sound that way.

You might try using ogr tools like "ogrinfo" to probe your original shapefiles
to see what the extents are before they're imported into grass. The
DNR data might not be in UTM coordinates or something, even if the rest of your
data is.

By the way, in an earlier email you said this this:

In that you were saying that
the use of -o allows v.in.ogr to assume that the projection is ok and
forces the points to plot within the projection's dimensions?

v.in.ogr -o doesn't force anything about plotting, per se. By default v.in.ogr
checks to see if the input data is in the same coordinate
system/datum/projection as the current location, and refuses to import the
ogr data source into a grass vector if it isn't. Your data has no
projection/coordinate system stored with it, so v.in.ogr assumes it's in the
wrong coordinates. But you *know* that the coordinates in your shapefile are
indeed in Zone 15 UTM coordinates, NAD83 datum, so you use "-o" to tell
v.in.ogr to skip the check and just import the data. No forcing involved ---
you're just telling v.in.ogr "it's OK, I know the coordinates in my input
files are in the right system, just read them, trust me." If, in fact, your
file is full of coordinates that are NOT in the right system and you use
v.in.ogr with "-o", well, then, you're on your own :). It sounds like you've
got a collection of data some of which is in UTM zone 15 NAD83, and some of
which might not be.

--
Tom Russo KM5VY SAR502 DM64ux http://www.swcp.com/~russo/
Tijeras, NM QRPL#1592 K2#398 SOC#236 AHTB#1
"The only thing you can do easily is be wrong, and that's hardly
  worth the effort." -- Norton Juster

7

On Mon, Jan 30, 2006 at 12:12:23PM -0600, we recorded a bogon-computron collision of the <kwythers@umn.edu> flavor, containing:

On Jan 30, 2006, at 10:05 AM, Tom Russo wrote:

>

Now here is what I think is going on with these damn files... tell me
if this makes sense to you. The coordinates in the .shp file from the
state are correct. This is why if I read the files in with v.in.ogr,
they can be displayed with d.vect on top of an outline of the state
is all looks fine. However, the y coordinates in the attribute file
(which come from the dbf file that came with the .shp file from the
state), are all messed up by some 4,700,206, plus or minus 0.5 (for
god knows what reason).

Ah! I understand now. The *attributes* of the points in the data file
have x and y columns, and those are weird. *sigh*

Here is what I'm wondering... if the .shp file has the correct
geography coordinates, is there some command that out write out those
coordinates? If so, I could substitute the good y coordinates for the
bad y coordinates in the attribute table. I assume if I use
v.out.ogr, that just writes the attributes from the attribute table
(including the weird y coordinates). true?

Yes, v.in.ogr is just taking the .dbf file's columns and importing them
into your postgresql database.

I am not really familiar with all the tricks in GRASS's db bag of tricks.
There is almost certainly some way to do it using db.execute, but you'll
have to get help from one of the folks on the GRASSLIST that actually knows
that stuff cold. It seems, though, that there must be a way for you to do
what you're trying to do using the non-geometry fields of the DB and the
good geometry data from the vectors, wihtout monkeying with the X and Y columns
of the attribute data.

--
Tom Russo KM5VY SAR502 DM64ux http://www.swcp.com/~russo/
Tijeras, NM QRPL#1592 K2#398 SOC#236 AHTB#1
"The only thing you can do easily is be wrong, and that's hardly
  worth the effort." -- Norton Juster