[GRASSLIST:5777] GRASS Novice Needs Help

John_Ronan · February 15, 2005, 10:48am

Morning,

I recently came across Grass when I posted on another list asking the question on how to convert scanned images into georeferenced images that could be used in xastir (www.xastir.org). So in fact Im not quite knowledgeable enough to call myself a novice even.

Thus far I've managed to convert the tiffs to unprojected lat/long rasters that xastir (www.xastir.org) can use and in fairness they seem to be pretty accurate. This was only possible by the extra-ordinary patience of Gerry Creager and Tom Russo. Who pretty much walked me through it step by step, thanks guys.

Anyways my problem now is trying to get rid of the 'fluff' that gets created around the edges by gdalwarp when it does the 'warping'. I've been attempting to use r.patch followed by r.proj to reproject the patched raster into another lat/long database but I've thus far not been successful.

I'm currently getting the error 'Input map is outside current region', here are the regions,

Firstly the Irish Grid.
projection: 99 (Transverse Mercator)
zone: 0
datum: ire65
ellipsoid: modif_airy
north: 52.3
south: 52.15
west: -7.4
east: -6.85
nsres: 0.15
ewres: 0.55
rows: 1
cols: 1

The Lat/Long region (bigger than the one above I think)
g.region -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 54N
south: 52N
west: 9W
east: 5W
nsres: 0:00:00.187198
ewres: 0:00:00.1872
rows: 38462
cols: 76923

So I've done a r.patch of two 'tiles' in the Irish Grid Database called ubermap. I exit grass, open it in the 'latlong' location and try to do r.proj

r.proj input=ubermap location=Ireland mapset=osmaps output=waterford1 method=nearest

Input Projection Parameters: +proj=tmerc +lat_0=53.5000000000 +lon_0=-8.0000000000 +k_0=1.0000350000 +x_0=200000.0000000000 +y_0=250000.0000000000 +a=6377340.189 +rf=299.3249646 +no_defs +towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15
Input Unit Factor: 1

Output Projection Parameters: +proj=latlong +a=6378137 +rf=298.257223563 +no_defs +towgs84=0.000,0.000,0.000
Output Unit Factor: 1
ERROR: Input map is outside current region

And I'm stuck... as I've really got no idea what I'm doing here any help would be appreciated.

Regards
de John
EI7IG

Ian_MacMillan1 · February 15, 2005, 6:49pm

John,
I think your problem is that your Irish grid is in transverse mercator, which means the units of measurement is meters, not degrees. Typically your coordinates will be measured in several hundred thousand to millions of meters. Latlong locations measure everything in degrees (like you have). Somehow, you will need to get the correct coordinates in meters for your irish grid first, and then r.proj into your latlong grid. Or mebbe you can try to just go straight into your latlong grid? Not sure if this is super helpful but mebbe this can help get you started.

G'luck
-Ian

On Feb 15, 2005, at 2:48 AM, John Ronan wrote:

Morning,

I recently came across Grass when I posted on another list asking the question on how to convert scanned images into georeferenced images that could be used in xastir (www.xastir.org). So in fact Im not quite knowledgeable enough to call myself a novice even.

Thus far I've managed to convert the tiffs to unprojected lat/long rasters that xastir (www.xastir.org) can use and in fairness they seem to be pretty accurate. This was only possible by the extra-ordinary patience of Gerry Creager and Tom Russo. Who pretty much walked me through it step by step, thanks guys.

Anyways my problem now is trying to get rid of the 'fluff' that gets created around the edges by gdalwarp when it does the 'warping'. I've been attempting to use r.patch followed by r.proj to reproject the patched raster into another lat/long database but I've thus far not been successful.

I'm currently getting the error 'Input map is outside current region', here are the regions,

Firstly the Irish Grid.
projection: 99 (Transverse Mercator)
zone: 0
datum: ire65
ellipsoid: modif_airy
north: 52.3
south: 52.15
west: -7.4
east: -6.85
nsres: 0.15
ewres: 0.55
rows: 1
cols: 1

The Lat/Long region (bigger than the one above I think)
g.region -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 54N
south: 52N
west: 9W
east: 5W
nsres: 0:00:00.187198
ewres: 0:00:00.1872
rows: 38462
cols: 76923

So I've done a r.patch of two 'tiles' in the Irish Grid Database called ubermap. I exit grass, open it in the 'latlong' location and try to do r.proj

r.proj input=ubermap location=Ireland mapset=osmaps output=waterford1 method=nearest

Input Projection Parameters: +proj=tmerc +lat_0=53.5000000000 +lon_0=-8.0000000000 +k_0=1.0000350000 +x_0=200000.0000000000 +y_0=250000.0000000000 +a=6377340.189 +rf=299.3249646 +no_defs +towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15
Input Unit Factor: 1

Output Projection Parameters: +proj=latlong +a=6378137 +rf=298.257223563 +no_defs +towgs84=0.000,0.000,0.000
Output Unit Factor: 1
ERROR: Input map is outside current region

And I'm stuck... as I've really got no idea what I'm doing here any help would be appreciated.

Regards
de John
EI7IG

John_Ronan · February 15, 2005, 8:23pm

On 15 Feb 2005, at 18:49, Ian MacMillan wrote:

John,
I think your problem is that your Irish grid is in transverse mercator, which means the units of measurement is meters, not degrees. Typically your coordinates will be measured in several hundred thousand to millions of meters. Latlong locations measure everything in degrees (like you have). Somehow, you will need to get the correct coordinates in meters for your irish grid first, and then r.proj into your latlong grid. Or mebbe you can try to just go straight into your latlong grid? Not sure if this is super helpful but mebbe this can help get you started.

Ok, I'm not following at all here...

Though I think I know what your getting at.. its not as simple converting from one projection system to another.

I've got the tiles in Irish Grid, and I can export them and gdalwarp them to unprojected lat/long which work fine (in xastir). As I can figure out the correct co-ordinates in Irish Grid (I had to in order to reproject the tiles in the first place) I'm not sure exactly how that helps me?

Of course the question remains, is this the correct way to 'patch' two 'tiles' together, which is my ultimate goal.

Regards
John
--
John Ronan <jronan at tssg dot org>, +353-51-302938
Telecommunications Software & Systems Group, http://www.tssg.org

Ian_MacMillan1 · February 15, 2005, 8:52pm

John, I will try to explain my reasoning a little more clearly. Now then, I could be wrong, but hopefully if someone else sees that I am leading you astray, they will chime in to this thread.

1) OK, I was trying to say that it sounds like your reprojection with r.proj is not going to work because your regions don't match in between your two locations. To reproject something from location 1 to location 2, you need to have the region in location 2 contain the coordinates of your raster from location 1. For example, you are trying to project a raster with bounding coord's of
north: 52.3
south: 52.15
west: -7.4
east: -6.85

into a location with coord's of

north: 54N
south: 52N
west: 9W
east: 5W

Normally this would work. However your coord's from the first location are not 52.3 degrees north to 52.15 degrees south, etc. They are 52.3 METERS north of the equator (or whatever the origin is for Irish grid), and -7.4 METERS west of the central meridian. Therefore your region in your latlong location does not contain the coordinates of your raster from your irish grid location. Following along these lines, I think your raster is georeferenced incorrectly in your Irish grid location.

Now then, if I understand what you are doing (taking a non-georeferenced raster image, and converting it into a georeferenced raster image such as a geotiff), here is the approach I would take.

Make an XY location, and import your nongeoreferenced rasters. Make a new location with projection parameters that match your final desired product. Go to the XY location, and use i.group, i.target, i.points, i.rectify. Open up your second location, and export your rasters with r.out.gdal.

2) It sounds like you are also having problems with the edges of your rasters, and you would like to crop them, and or patch them together to get rid of these blemishes. If you use r.patch, make sure that the first map listed in input= is the top map. If you want to crop the maps, just set your region to the raster with g.region rast=your_rast. Then run d.zoom, and zoom into the region that you want to be the new 'cropped' raster. Exit d.zoom, and run r.mapcalc your_newrast=your_rast. This will give you a cropped map based on the extent of the region settings.

Hope this helps
-Ian

On Feb 15, 2005, at 12:23 PM, John Ronan wrote:

On 15 Feb 2005, at 18:49, Ian MacMillan wrote:

John,
I think your problem is that your Irish grid is in transverse mercator, which means the units of measurement is meters, not degrees. Typically your coordinates will be measured in several hundred thousand to millions of meters. Latlong locations measure everything in degrees (like you have). Somehow, you will need to get the correct coordinates in meters for your irish grid first, and then r.proj into your latlong grid. Or mebbe you can try to just go straight into your latlong grid? Not sure if this is super helpful but mebbe this can help get you started.

Ok, I'm not following at all here...

Though I think I know what your getting at.. its not as simple converting from one projection system to another.

I've got the tiles in Irish Grid, and I can export them and gdalwarp them to unprojected lat/long which work fine (in xastir). As I can figure out the correct co-ordinates in Irish Grid (I had to in order to reproject the tiles in the first place) I'm not sure exactly how that helps me?

Of course the question remains, is this the correct way to 'patch' two 'tiles' together, which is my ultimate goal.

Regards
John
--
John Ronan <jronan at tssg dot org>, +353-51-302938
Telecommunications Software & Systems Group, http://www.tssg.org

John_Ronan · February 15, 2005, 9:22pm

On 15 Feb 2005, at 20:52, Ian MacMillan wrote:

John, I will try to explain my reasoning a little more clearly. Now then, I could be wrong, but hopefully if someone else sees that I am leading you astray, they will chime in to this thread.

1) OK, I was trying to say that it sounds like your reprojection with r.proj is not going to work because your regions don't match in between your two locations. To reproject something from location 1 to location 2, you need to have the region in location 2 contain the coordinates of your raster from location 1. For example, you are trying to project a raster with bounding coord's of
north: 52.3
south: 52.15
west: -7.4
east: -6.85

into a location with coord's of

north: 54N
south: 52N
west: 9W
east: 5W

Normally this would work. However your coord's from the first location are not 52.3 degrees north to 52.15 degrees south, etc. They are 52.3 METERS north of the equator (or whatever the origin is for Irish grid), and -7.4 METERS west of the central meridian. Therefore your region in your latlong location does not contain the
coordinates of your raster from your irish grid location. Following along these lines, I think your raster is georeferenced incorrectly in your Irish grid location.

Doh!

Sorry about that.

Now then, if I understand what you are doing (taking a non-georeferenced raster image, and converting it into a georeferenced raster image such as a geotiff), here is the approach I would take.

Make an XY location, and import your nongeoreferenced rasters. Make a new location with projection parameters that match your final desired product. Go to the XY location, and use i.group, i.target, i.points, i.rectify. Open up your second location, and export your rasters with r.out.gdal.

Thats what I've done... and they seem to be fine.

2) It sounds like you are also having problems with the edges of your rasters, and you would like to crop them, and or patch them together to get rid of these blemishes. If you use r.patch, make sure that the first map listed in input= is the top map. If you want to crop the maps, just set your region to the raster with g.region rast=your_rast. Then run d.zoom, and zoom into the region that you want to be the new 'cropped' raster. Exit d.zoom, and run r.mapcalc your_newrast=your_rast. This will give you a cropped map based on the extent of the region settings.

Ok this is the bit I'm having trouble with. I don't want to crop(yet), its just where they meet I'm getting rubbish appearing from the edge of one tile over the other as gdalwarp (I assume) doesn't have any data to work with to fill in the missing bits so I just get a white bit where the two tiles meet. This is what I'm trying to remove. As both tiles seem to be going right to the edge of one another (when I zoom in it looks like GRASS hasn't interpolated and I can see 'through' one tile to the other)

If I can get it working with two tiles first, I'll progress from there

Hope this helps

Yes it has indeed. I'll remember the cropping bit as I may need to do that in the future.

Regards
John

--
John Ronan <jronan at tssg dot org>, +353-51-302938
Telecommunications Software & Systems Group, http://www.tssg.org

H_B · February 16, 2005, 1:05am

I'm currently getting the error 'Input map is outside current region',

here are the regions,

Firstly the Irish Grid.
projection: 99 (Transverse Mercator)
zone: 0
datum: ire65
ellipsoid: modif_airy
north: 52.3
south: 52.15
west: -7.4
east: -6.85

[as noted, there are wrong; what does "gdalinfo" say about the TIFF?
are you sure the TIFF isn't starting out in lat/lon already?]

nsres: 0.15
ewres: 0.55
rows: 1
cols: 1

so resolution is 15cm x 55cm and everything is contained in one single
data cell.

The Lat/Long region (bigger than the one above I think)
g.region -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 54N
south: 52N
west: 9W
east: 5W
nsres: 0:00:00.187198
ewres: 0:00:00.1872
rows: 38462
cols: 76923

38462 rows X 76923 columns! You've chosen a very large area of land to
map out a 5m resolution!

To patch:

g.region rast=map1,map2
## strech region to contain both maps

r.patch in=map1,map2 out=map1_2

And I'm stuck... as I've really got no idea what I'm doing here any
help would be appreciated.

Have you tried the tutorials?
http://grass.ibiblio.org/ and follow the links... also have a look on
the Wiki.

Hamish

John_Ronan · February 16, 2005, 10:01am

[as noted, there are wrong; what does "gdalinfo" say about the TIFF?
are you sure the TIFF isn't starting out in lat/lon already?]

Ok, the tiff generated by r.out.gdal give

gdalinfo os2410poly1.tif
Driver: GTiff/GeoTIFF
Size is 4724, 4724
Coordinate System is `'
Origin = (-7.415917,52.330316)
Pixel Size = (0.000063,-0.000038)
Corner Coordinates:
Upper Left ( -7.4159169, 52.3303161)
Lower Left ( -7.4159169, 52.1487398)
Upper Right ( -7.1200431, 52.3303161)
Lower Right ( -7.1200431, 52.1487398)
Center ( -7.2679800, 52.2395280)

nsres: 0.15
ewres: 0.55
rows: 1
cols: 1

so resolution is 15cm x 55cm and everything is contained in one single
data cell.

When you explain it like that, it doesn't make much sense does it? Hmmm

[snip]

38462 rows X 76923 columns! You've chosen a very large area of land to
map out a 5m resolution!

Well I was trying to make it large enough

To patch:

g.region rast=map1,map2
## strech region to contain both maps

Ok, something funny going on here thats the output from g.region -p from both maps individually

GRASS 5.7.0:~ > g.region -p
projection: 99 (Transverse Mercator)
zone: 0
datum: ire65
ellipsoid: modif_airy
north: 52.33030351
south: 52.14872735
west: -7.41558677
east: -7.11971255
nsres: 0.00003844
ewres: 0.00006263
rows: 4724
cols: 4724

GRASS 5.7.0:~ > g.region -p
projection: 99 (Transverse Mercator)
zone: 0
datum: ire65
ellipsoid: modif_airy
north: 52.32853668
south: 52.14617722
west: -7.12318884
east: -6.8268459
nsres: 0.0000386
ewres: 0.00006273
rows: 4724
cols: 4724

Which doesn't at all agree with Tom's follow on email, Hmmm.

I suspect I've done something else stupid. When GRASS asked me for the tie points for the rasters, I supplied the eastings northings in latitude and longitude. Which I'm assuming now given that my co-ordinate system is in Meters that this is completely wrong. That wouldn't explain how they seem to display correctly though.

Tom, how did you work out the tie points in meters?

Regards
John
--
John Ronan <jronan at tssg dot org>, +353-51-302938
Telecommunications Software & Systems Group, http://www.tssg.org

Tom_Russo · February 16, 2005, 4:45pm

On Wed, Feb 16, 2005 at 10:01:54AM +0000, we recorded a bogon-computron collision of the <jronan@tssg.org> flavor, containing:

Ok, something funny going on here thats the output from g.region -p
from both maps individually

GRASS 5.7.0:~ > g.region -p
projection: 99 (Transverse Mercator)
zone: 0
datum: ire65
ellipsoid: modif_airy
north: 52.33030351
south: 52.14872735
west: -7.41558677
east: -7.11971255
nsres: 0.00003844
ewres: 0.00006263
rows: 4724
cols: 4724

[...]

Which doesn't at all agree with Tom's follow on email, Hmmm.

I suspect I've done something else stupid. When GRASS asked me for the
tie points for the rasters, I supplied the eastings northings in
latitude and longitude. Which I'm assuming now given that my
co-ordinate system is in Meters that this is completely wrong. That
wouldn't explain how they seem to display correctly though.

Yes, you need to give your tie-points in the northings and eastings of the
Irish Grid System which are in meters.

Tom, how did you work out the tie points in meters?

Your maps have a grid that show the northings and eastings --- but obfuscated
a little. Just as the US Military abbreviates UTMs by turning the first
few digits of each coordinate into a 2-letter grid square, and compresses
the remaining digits into a single 6-digit code, the Irish Grid system
abbreviates the 6-digit northings and eastings into a single letter and two
3-digit partial coordinates. Then on maps they throw the least significant
digit away when labeling the grid lines.

Your OSi maps have a 1000m grid with intersections marked as, for example,
"41 11"
which is an abbreviation of the coordinates
"S 410 110"
(you told me that) as recognized by the Irish Grid calculator you first
pointed me at. With a little digging I determined that the "S" was an
encoding of the 100Km digit of the northings and eastings, and that "410" and
"110" are the 10k, 1k, and 100m digits of the coordinates.

In the case of the grid intersection 41/11 or abbreviated coordinate
"S 410 110" the correct coordinates in meters are 241000, 111000.

For other lettered grids you'll have to work out what the first digits are
supposed to be, then you can map the 2-digit grid intersctions into 6 digit
eastings and northings quite easily. It is these coordinates you need to use
to label tie points in i.points.

To convert between known Irish Grid coordinates (in meters) to lat/lon
(for sanity-checking purposes), you can use the cs2cs command:

    cs2cs +proj=tmerc +lat_0=53.5000000000 +lon_0=-8.0000000000 \
        +k_0=1.0000350000 +x_0=200000.0000000000 +y_0=250000.0000000000 \
        +a=6377340.189 +rf=299.3249646 +no_defs \
        +towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15 \
    +to +proj=latlong +datum=WGS84

To do it the other way, use:
    cs2cs +proj=latlong +datum=WGS84 \
     +to +proj=tmerc +lat_0=53.5000000000 +lon_0=-8.0000000000 \
        +k_0=1.0000350000 +x_0=200000.0000000000 +y_0=250000.0000000000 \
        +a=6377340.189 +rf=299.3249646 +no_defs \
        +towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15 \

This is how I managed to get the fact that "S" meant "put a 2 in front of
the easting and a 1 in front of the northing" --- I used the converter web
site you told me about to convert S 410 110 into lat/lon, then used
cs2cs to convert it back into meters in Irish Grid. When I say "sanity
checking" I mean double-checking in the final quality control step, e.g. plug
"241000 111000" into the first cs2cs command and you should find the
lat/lon that puts you at the 41/11 grid intersection on the waterford/tramore
map I produced back in October.

Hope that helps. Glad you posted here, because I was assuming you were
using the northings and eastings in Irish Grid when I was trying to figure
out your other questions.

--
Tom Russo KM5VY SAR502 DM64ux http://www.swcp.com/~russo/
Tijeras, NM QRPL#1592 K2#398 SOC#236 AHTB#1
"When life gives you lemons, find someone with a paper cut."