[GRASS-user] trying to get unprojected raster data into a projected location

GRASS GRASS user
,
1

I received a set of data that was developed in ArgGIS, including shapefiles and a .sid raster file. The shapefiles came with projection information, but the .sid file did not (and it extends out further than the vector info on the edges). So I created an x,y location for the .sid file and imported it using GDAL. It imported fine; the satphoto looks good. I also created a UTM location for the shapefiles using their proj info, and then extended the region and increased the resolution based on the extents of the raster (all units are in meters).
Now I want to reproject the raster from the x,y location into the UTM location. Both are in meters, both have the same N,S,W,E extents, both have the same resolution.

I cannot figure out the steps, based on the documentation I have found through 3 hours of searching.
First, I tried to create a new long-lat location using the same extents, but in GRASS, creating a location manually seems to default to d-m-s coordinates and does not allow for meters. So I cannot input my extents, false northing, false easting, and resolution.
Second, I tried using g.setproj within the x,y location but that is not allowed.

So I cannot seem to generate projection information for this raster that I have imported, and therefore I cannot reproject it into the UTM location with the same extents, resolution, and units.

It seems to me that this problem should be common, and the solution should be elementary. It should be easy to find documentation on how to resolve this problem in introductory tutorials. I can imagine that many GRASS users receive ESRI-generated data with missing or partial projection information. Please advise.

Pietro Calogero
Mobile: +93 (0) 799 02 00 27

http://www.calogero.us/Kabul/2007/

Skype: pietrocalogero

2

Pietro Calogero wrote:

I received a set of data that was developed in ArgGIS, including
shapefiles and a .sid raster file. The shapefiles came with
projection information, but the .sid file did not (and it extends
out further than the vector info on the edges)

This shouldn't be a big problem. If you are sure your SID raster is in
the same projection that your shapefiles just use flag -o in r.in.gdal
which will override the projection check.

Let me know if that helps.

Maciek

3

Pietro,

See below

On 9/20/07 8:56 PM, "Pietro Calogero" <pietro@calogero.us> wrote:

I received a set of data that was developed in ArgGIS, including shapefiles
and a .sid raster file. The shapefiles came with projection information, but
the .sid file did not (and it extends out further than the vector info on the
edges). So I created an x,y location for the .sid file and imported it using
GDAL. It imported fine; the satphoto looks good. I also created a UTM location
for the shapefiles using their proj info, and then extended the region and
increased the resolution based on the extents of the raster (all units are in
meters).

If you KNOW the projection of the .sid file AND it IS projected (just
doesn't have a header with projection information) you can initially create
a proper location and simply import the file overriding its (lack of)
projection information. This is a check box in the GUI and the -o flag on
the command line. No need to reproject the file.

Now I want to reproject the raster from the x,y location into the UTM
location. Both are in meters, both have the same N,S,W,E extents, both have
the same resolution.

If you have the file in an xy location and want project it into a UTM
location, the simplest way is to use the georectification module in the GUI.

Start GRASS in the proper UTM location. Make sure that you set the region
extents and resolution to match the image extents and native resolution for
best results. Use g.region to do this.

First, display one of your vector maps in the normal way.

Start the georectification module under the file menu.

1) select the location/mapset of the file you want to georectify
2) create a georectifcation "group" if you have not already defined one.
Your group will contain only 1 file
3) select that group
4) select a map/image to display for georectification. In your case, it is
your single image.
5) start georectification

If you don't know the correspondence between xy points and desired UTM
equivalents, you could interactively click a point on your image and a
corresponding point on your vector to get that information for a series of
points.

However, in your case, it sounds like you already know that there is a 1:1
correspondence between x/y values and UTM east/north and your map is already
rectified (i.e., doesn't need any kind of warping). If so, you can just type
the corner coordinates into the xy box (separate by commas as in
xcoord,ycoord) and the same number into the eastnorth box (format as
east,north). Do this for the four corners, pick 1st order rectification and
click the georectify button.

Hope this helps.

Michael

__________________________________________
Michael Barton, Professor of Anthropology
Director of Graduate Studies
School of Human Evolution & Social Change
Center for Social Dynamics & Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

4

Thank you Michael. Good to know about that capability of the georectify tool.
Yesterday Maciej emailed me a different solution, which also worked. I thought my question was so elementary that users on this list would not want it posted. But he suggested that I do that, for the benefit of other newbies, I suppose. So here is a summa:

  1. The shapefiles had .proj files associated with them. I used one to create a UTM location in GRASS.
  2. The .sid raster had no projection info attached to it (so far as I can tell). I used it to create an x,y location.
  3. The extents of the x.y location were broader, and the cell resolution much higher than the UTM location. I used g.region to increase the extents and resolution of the UTM projection to match the x,y location (and presumably the .sid raster).
  4. I imported the .sid raster directly into the UTM location. Actually, I only imported the red band (band=1) because that alone was 1.7 GB and ground features resolved most clearly in that particular band.
  5. The imported vectors and raster line up, so at least I have replicated the registration of the original ArcGIS coverage.

So in this procedure I only created the x,y location to extract information from the .sid file: its extents and resolution. I never pulled data into the x,y location. Thank God, because that raster takes 2 hours to import into any location!

Now struggling with v.overlay of intersecting buffers over polygons…grrrr…
Ciao,

Pietro Calogero

http://www.calogero.us/Kabul/2007/

On Sep 22, 2007, at 9:00 PM, Michael Barton wrote:
Pietro, see below

On 9/20/07 8:56 PM, “Pietro Calogero” <pietro@calogero.us> wrote:

I received a set of data that was developed in ArgGIS, including shapefiles
and a .sid raster file. The shapefiles came with projection information, but
the .sid file did not (and it extends out further than the vector info on the
edges). So I created an x,y location for the .sid file and imported it using
GDAL. It imported fine; the satphoto looks good. I also created a UTM location
for the shapefiles using their proj info, and then extended the region and
increased the resolution based on the extents of the raster (all units are in
meters).

If you KNOW the projection of the .sid file AND it IS projected (just
doesn’t have a header with projection information) you can initially create
a proper location and simply import the file overriding its (lack of)
projection information. This is a check box in the GUI and the -o flag on
the command line. No need to reproject the file.

Now I want to reproject the raster from the x,y location into the UTM
location. Both are in meters, both have the same N,S,W,E extents, both have
the same resolution.

If you have the file in an xy location and want project it into a UTM
location, the simplest way is to use the georectification module in the GUI.

Start GRASS in the proper UTM location. Make sure that you set the region
extents and resolution to match the image extents and native resolution for
best results. Use g.region to do this.

First, display one of your vector maps in the normal way.

Start the georectification module under the file menu.

  1. select the location/mapset of the file you want to georectify
  2. create a georectifcation “group” if you have not already defined one.
    Your group will contain only 1 file
  3. select that group
  4. select a map/image to display for georectification. In your case, it is
    your single image.
  5. start georectification

If you don’t know the correspondence between xy points and desired UTM
equivalents, you could interactively click a point on your image and a
corresponding point on your vector to get that information for a series of
points.

However, in your case, it sounds like you already know that there is a 1:1
correspondence between x/y values and UTM east/north and your map is already
rectified (i.e., doesn’t need any kind of warping). If so, you can just type
the corner coordinates into the xy box (separate by commas as in
xcoord,ycoord) and the same number into the eastnorth box (format as
east,north). Do this for the four corners, pick 1st order rectification and
click the georectify button.

Hope this helps.

Michael

Michael Barton, Professor of Anthropology
Director of Graduate Studies
School of Human Evolution & Social Change
Center for Social Dynamics & Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

5

Hi,

You could use gdalinfo utility to get info about raster resolution and
size thus importing into GRASS X,Y location just to see it's
resolution is unnecessary.

Maris.

2007/9/22, Pietro Calogero <pietro@calogero.us>:

So in this procedure I only created the x,y location to extract
information from the .sid file: its extents and resolution. I never
pulled data into the x,y location. Thank God, because that raster
takes 2 hours to import into any location!

6

5. The imported vectors and raster line up, so at least I have
replicated the registration of the original ArcGIS coverage.

'r.in.gdal -o' doesn't pay any attention to the g.region settings, it just uses
the data in the file.

If it imported as south and west = 0, and north and east = number of pixels,
with a resolution of 1, r.region can be used to reset that if the correct
values are known. (that won't work in lat/lon location where the number of
vertical pixels exceeds 90 (ie can't be north of the north pole), then the
import fails [devs: in that case should we rescale the import in r.in.gdal to
y=0-1?])

So in this procedure I only created the x,y location to extract
information from the .sid file: its extents and resolution. I never
pulled data into the x,y location. Thank God, because that raster
takes 2 hours to import into any location!

gdalinfo should be able to give you that info without importing it.
gdal_translate might let you assign georeferencing info to the image if you
wanted to do that, but as a new file.

Hamish

____________________________________________________________________________________
Pinpoint customers who are looking for what you sell.
http://searchmarketing.yahoo.com/

7

Thank you. GDAL has capabilities I was not aware of. Part of the challenge here is that some ArcView users do not generate .proj files even with shapefiles (not automatic in 3.2), and generally the GIS users around here operate in metric units, so using the Lat-Lon location cannot work. And I have not figured out how to change an x,y location so that data can be reprojected from it into UTM.
In a case where I do not have a match as described earlier, I guess I will need to use the georeference tool to generate proj data before importing. I have done that easily in QGis with 6 MB Google Earth mashup TIFs, so I suppose the method is similar for unreferenced shapefiles. But I’m not sure how QGIS handles a 5GB .sid file.

Pietro Calogero

http://www.calogero.us/Kabul/2007/

On Sep 23, 2007, at 8:37 AM, Hamish wrote:

  1. The imported vectors and raster line up, so at least I have
    replicated the registration of the original ArcGIS coverage.

‘r.in.gdal -o’ doesn’t pay any attention to the g.region settings, it just uses
the data in the file.

If it imported as south and west = 0, and north and east = number of pixels,
with a resolution of 1, r.region can be used to reset that if the correct
values are known. (that won’t work in lat/lon location where the number of
vertical pixels exceeds 90 (ie can’t be north of the north pole), then the
import fails [devs: in that case should we rescale the import in r.in.gdal to
y=0-1?])

So in this procedure I only created the x,y location to extract
information from the .sid file: its extents and resolution. I never
pulled data into the x,y location. Thank God, because that raster
takes 2 hours to import into any location!

gdalinfo should be able to give you that info without importing it.
gdal_translate might let you assign georeferencing info to the image if you
wanted to do that, but as a new file.

Hamish

Pinpoint customers who are looking for what you sell.
http://searchmarketing.yahoo.com/