[GRASS-user] Help: Converting a raster map between locations (from wgs84 lat / long to UTM)

Dear friends,

I am trying to convert a GTOPO30 DEM (raster) map (for UK) from lat / long wgs
84 to UTM for UK.

I am running into rserious problems. I am using the procedure as described in
the book open source gis a grass gis approach"

1) in the WGS84 lat long location (ukll) I resize the map using d.zoom
2) I save the new region as well as the new sub-raster (for scotland)
3) I open grass in the UTM location
4) I try to import using

r.proj in=gtopo30_cs out=gtopo30_cs_utm location=ukll mapset=PERMANENT
method=cubic --o

All I can get is: ERROR: Input map is outside current region

What am I doing wrong? I have tried so many different alternative solution.

Best,
--
Corrado Topi

Global Climate Change & Biodiversity Indicators
Area 18,Department of Biology
University of York, York, YO10 5YW, UK
Phone: + 44 (0) 1904 328645, E-mail: ct529@york.ac.uk

It's because your region in the target location is wrong.

Here is a little trick:

In your lat/lon location, after you set your region, run the command
v.in.region (check correct syntax). This will create a bounding box of
the region.

Then, go to your target location (UTM) and use the command v.proj to
import the vector you just created.

Now, set your region to the imported vector with g.region
vect=name_of_vector_here

Finally, run r.proj!!

The only problem at which I don't know any exact method is setting the
target resolution correctly... I know it's done with g.region res=XX
but, the conversion from lat/long to meters is not a straight one...

Cheers
Daniel

On Wed, Oct 22, 2008 at 2:44 PM, Corrado <ct529@york.ac.uk> wrote:

Dear friends,

I am trying to convert a GTOPO30 DEM (raster) map (for UK) from lat / long wgs
84 to UTM for UK.

I am running into rserious problems. I am using the procedure as described in
the book open source gis a grass gis approach"

1) in the WGS84 lat long location (ukll) I resize the map using d.zoom
2) I save the new region as well as the new sub-raster (for scotland)
3) I open grass in the UTM location
4) I try to import using

r.proj in=gtopo30_cs out=gtopo30_cs_utm location=ukll mapset=PERMANENT
method=cubic --o

All I can get is: ERROR: Input map is outside current region

What am I doing wrong? I have tried so many different alternative solution.

Best,
--
Corrado Topi

Global Climate Change & Biodiversity Indicators
Area 18,Department of Biology
University of York, York, YO10 5YW, UK
Phone: + 44 (0) 1904 328645, E-mail: ct529@york.ac.uk
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

On Wed, 2008-10-22 at 17:44 +0100, Corrado wrote:

Dear friends,

I am trying to convert a GTOPO30 DEM (raster) map (for UK) from lat / long wgs
84 to UTM for UK.

I am running into rserious problems. I am using the procedure as described in
the book open source gis a grass gis approach"

1) in the WGS84 lat long location (ukll) I resize the map using d.zoom
2) I save the new region as well as the new sub-raster (for scotland)
3) I open grass in the UTM location
4) I try to import using

r.proj in=gtopo30_cs out=gtopo30_cs_utm location=ukll mapset=PERMANENT
method=cubic --o

All I can get is: ERROR: Input map is outside current region

What am I doing wrong? I have tried so many different alternative solution.

Best,

Hi Corrado.

You need to check the extent of the map you want to reproject in the
initial location/mapset (check with g.region -p) and then set those (or
bigger?) in your target location/mapset (with g.region e= s= w= n= ).

Then re-"r.proj" and it should go.

Cheers, Nikos

P.S. If you are interested why it's not "automatic" you can read more in
this "trac ticket" http://trac.osgeo.org/grass/ticket/328

On Wed, 2008-10-22 at 15:24 -0200, Daniel Victoria wrote:

It's because your region in the target location is wrong.

Here is a little trick:

In your lat/lon location, after you set your region, run the command
v.in.region (check correct syntax). This will create a bounding box of
the region.

Then, go to your target location (UTM) and use the command v.proj to
import the vector you just created.

Now, set your region to the imported vector with g.region
vect=name_of_vector_here

Finally, run r.proj!!

The only problem at which I don't know any exact method is setting the
target resolution correctly... I know it's done with g.region res=XX
but, the conversion from lat/long to meters is not a straight one...

Cheers
Daniel

On Wed, Oct 22, 2008 at 2:44 PM, Corrado <ct529@york.ac.uk> wrote:
> Dear friends,
>
> I am trying to convert a GTOPO30 DEM (raster) map (for UK) from lat / long wgs
> 84 to UTM for UK.
>
> I am running into rserious problems. I am using the procedure as described in
> the book open source gis a grass gis approach"
>
> 1) in the WGS84 lat long location (ukll) I resize the map using d.zoom
> 2) I save the new region as well as the new sub-raster (for scotland)
> 3) I open grass in the UTM location
> 4) I try to import using
>
> r.proj in=gtopo30_cs out=gtopo30_cs_utm location=ukll mapset=PERMANENT
> method=cubic --o
>
> All I can get is: ERROR: Input map is outside current region
>
> What am I doing wrong? I have tried so many different alternative solution.
>
> Best,
> --
> Corrado Topi

OK, maybe your way is better Daniel :slight_smile:

Greets, Nikos

Thanks Nikos, thanks Daniel!

Unfortunately, my map is a raster map, not a vector .... can I still use
v.in.region and v.proj?

Thanks

On Wednesday 22 October 2008 18:27:23 Nikos Alexandris wrote:

On Wed, 2008-10-22 at 15:24 -0200, Daniel Victoria wrote:
> It's because your region in the target location is wrong.
>
> Here is a little trick:
>
> In your lat/lon location, after you set your region, run the command
> v.in.region (check correct syntax). This will create a bounding box of
> the region.
>
> Then, go to your target location (UTM) and use the command v.proj to
> import the vector you just created.
>
> Now, set your region to the imported vector with g.region
> vect=name_of_vector_here
>
> Finally, run r.proj!!
>
> The only problem at which I don't know any exact method is setting the
> target resolution correctly... I know it's done with g.region res=XX
> but, the conversion from lat/long to meters is not a straight one...
>
> Cheers
> Daniel
>
> On Wed, Oct 22, 2008 at 2:44 PM, Corrado <ct529@york.ac.uk> wrote:
> > Dear friends,
> >
> > I am trying to convert a GTOPO30 DEM (raster) map (for UK) from lat /
> > long wgs 84 to UTM for UK.
> >
> > I am running into rserious problems. I am using the procedure as
> > described in the book open source gis a grass gis approach"
> >
> > 1) in the WGS84 lat long location (ukll) I resize the map using d.zoom
> > 2) I save the new region as well as the new sub-raster (for scotland)
> > 3) I open grass in the UTM location
> > 4) I try to import using
> >
> > r.proj in=gtopo30_cs out=gtopo30_cs_utm location=ukll mapset=PERMANENT
> > method=cubic --o
> >
> > All I can get is: ERROR: Input map is outside current region
> >
> > What am I doing wrong? I have tried so many different alternative
> > solution.
> >
> > Best,
> > --
> > Corrado Topi

OK, maybe your way is better Daniel :slight_smile:

Greets, Nikos

--
Corrado Topi

Global Climate Change & Biodiversity Indicators
Area 18,Department of Biology
University of York, York, YO10 5YW, UK
Phone: + 44 (0) 1904 328645, E-mail: ct529@york.ac.uk

On Thu, 2008-10-23 at 10:40 +0100, Corrado wrote:

Thanks Nikos, thanks Daniel!

Unfortunately, my map is a raster map, not a vector .... can I still
use
v.in.region and v.proj?

Thanks

[...]

Hi Corrado!

For sure!! You can use v.in.region anytime since it deals with the
current region and not with a specific map.

Regards, Nikos

Corrado,

What v.in.region does is create a vector that surrounds your current
region. Then, when you run v.proj, the vector covering the region is
brought to your target location so you know where your projected
raster should be at...

The "manual labor" after that is setting the resolution correctly. You
could then use some rule of thumb like 3arcsec ~= 90 m or, try to set
the number of columns/rows in your target region similar to your
lat/lon region. But this latter option does not work so well if your
raster is to be tilted a lot...

On an "additional" question: How do people set the resolution
correctly when transforming from lat/lon to meter and vice-versa? Are
there any other methods besides the ones I mentioned?

Cheers
Daniel

On Thu, Oct 23, 2008 at 8:16 AM, Nikos Alexandris
<nikos.alexandris@felis.uni-freiburg.de> wrote:

On Thu, 2008-10-23 at 10:40 +0100, Corrado wrote:

Thanks Nikos, thanks Daniel!

Unfortunately, my map is a raster map, not a vector .... can I still
use
v.in.region and v.proj?

Thanks

[...]

Hi Corrado!

For sure!! You can use v.in.region anytime since it deals with the
current region and not with a specific map.

Regards, Nikos

On Thursday 23 October 2008 11:38:47 Daniel Victoria wrote:
Thanks Nikos and Daniel again!

I was seriously hoping that GRASS would have managed the conversion for
me .... I do not think we could actually define the resolution by hand (it
would not really be accepted in a paper)!

Corrado,

What v.in.region does is create a vector that surrounds your current
region. Then, when you run v.proj, the vector covering the region is
brought to your target location so you know where your projected
raster should be at...

The "manual labor" after that is setting the resolution correctly. You
could then use some rule of thumb like 3arcsec ~= 90 m or, try to set
the number of columns/rows in your target region similar to your
lat/lon region. But this latter option does not work so well if your
raster is to be tilted a lot...

On an "additional" question: How do people set the resolution
correctly when transforming from lat/lon to meter and vice-versa? Are
there any other methods besides the ones I mentioned?

Cheers
Daniel

On Thu, Oct 23, 2008 at 8:16 AM, Nikos Alexandris

<nikos.alexandris@felis.uni-freiburg.de> wrote:
> On Thu, 2008-10-23 at 10:40 +0100, Corrado wrote:
>> Thanks Nikos, thanks Daniel!
>>
>> Unfortunately, my map is a raster map, not a vector .... can I still
>> use
>> v.in.region and v.proj?
>>
>> Thanks
>
> [...]
>
> Hi Corrado!
>
> For sure!! You can use v.in.region anytime since it deals with the
> current region and not with a specific map.
>
> Regards, Nikos

--
Corrado Topi

Global Climate Change & Biodiversity Indicators
Area 18,Department of Biology
University of York, York, YO10 5YW, UK
Phone: + 44 (0) 1904 328645, E-mail: ct529@york.ac.uk

Corrado,

I was seriously hoping that GRASS would have managed the conversion for me ....

Yes, I agree that's annoying

I do not think we could actually define the resolution by hand (it would not really be accepted in a paper)!

Why not?
http://www.uwgb.edu/dutchs/UsefulData/HowUseExcel.HTM

Marco

Dear Daniel, dear Nikos, dear list,

I have tried everything I could. I read the GRASS books chapters and the
documentation. It does not work.

The problem is that the files I have are in WGS84 unprojected. See here:

http://www.worldclim.org/format.htm

Whichever procedure I use, whther Daniel's or the book (for example: page 47),
the problem is always the same.

I create a location + mapset + map with the worldclim data (location01). I
load the original data, and I call the raster map "precipitation". I generate
the correct box with v.region and I call it "test". I create a new location +
mapset (utm + osgb), called location02. I run v.proj in location B pointing
at map "test"

v.proj input=test location=location01 mapset=PERMANENT

I have the same problem:

PROJ_INFO file not found for location location01

But of course the projection is not defined, the data are unprojected!

I am slowly getting desperate .... I have been working on it several days and
cannot get a sensible result. I have tried with different datum
transformation parameters for osgb36 and trying to project .... but the
results are very unreliable ....

Do you know what should I do?

On Thursday 23 October 2008 11:38:47 Daniel Victoria wrote:

Corrado,

What v.in.region does is create a vector that surrounds your current
region. Then, when you run v.proj, the vector covering the region is
brought to your target location so you know where your projected
raster should be at...

The "manual labor" after that is setting the resolution correctly. You
could then use some rule of thumb like 3arcsec ~= 90 m or, try to set
the number of columns/rows in your target region similar to your
lat/lon region. But this latter option does not work so well if your
raster is to be tilted a lot...

On an "additional" question: How do people set the resolution
correctly when transforming from lat/lon to meter and vice-versa? Are
there any other methods besides the ones I mentioned?

Cheers
Daniel

On Thu, Oct 23, 2008 at 8:16 AM, Nikos Alexandris

<nikos.alexandris@felis.uni-freiburg.de> wrote:
> On Thu, 2008-10-23 at 10:40 +0100, Corrado wrote:
>> Thanks Nikos, thanks Daniel!
>>
>> Unfortunately, my map is a raster map, not a vector .... can I still
>> use
>> v.in.region and v.proj?
>>
>> Thanks
>
> [...]
>
> Hi Corrado!
>
> For sure!! You can use v.in.region anytime since it deals with the
> current region and not with a specific map.
>
> Regards, Nikos

--
Corrado Topi

Global Climate Change & Biodiversity Indicators
Area 18,Department of Biology
University of York, York, YO10 5YW, UK
Phone: + 44 (0) 1904 328645, E-mail: ct529@york.ac.uk

Even though the projection is not defined, the data still has some
information it is WGS84 which has a datum/coordinate system so you do
need to define it, the EPSG code is 4326.

You can't reproject something from nothing, it needs at least a
coordinate system, otherwise it doesn't know what point on your current
map should be moved and to where it should be moved.

When you configured your first location you should have specified the
EPSG as 4326.

Does that get you further?
Alex

Corrado wrote:

Dear Daniel, dear Nikos, dear list,

I have tried everything I could. I read the GRASS books chapters and the
documentation. It does not work.

The problem is that the files I have are in WGS84 unprojected. See here:

http://www.worldclim.org/format.htm

Whichever procedure I use, whther Daniel's or the book (for example: page 47),
the problem is always the same.

I create a location + mapset + map with the worldclim data (location01). I
load the original data, and I call the raster map "precipitation". I generate
the correct box with v.region and I call it "test". I create a new location +
mapset (utm + osgb), called location02. I run v.proj in location B pointing
at map "test"

v.proj input=test location=location01 mapset=PERMANENT

I have the same problem:

PROJ_INFO file not found for location location01

But of course the projection is not defined, the data are unprojected!

I am slowly getting desperate .... I have been working on it several days and
cannot get a sensible result. I have tried with different datum
transformation parameters for osgb36 and trying to project .... but the
results are very unreliable ....

Do you know what should I do?

On Thursday 23 October 2008 11:38:47 Daniel Victoria wrote:

Corrado,

What v.in.region does is create a vector that surrounds your current
region. Then, when you run v.proj, the vector covering the region is
brought to your target location so you know where your projected
raster should be at...

The "manual labor" after that is setting the resolution correctly. You
could then use some rule of thumb like 3arcsec ~= 90 m or, try to set
the number of columns/rows in your target region similar to your
lat/lon region. But this latter option does not work so well if your
raster is to be tilted a lot...

On an "additional" question: How do people set the resolution
correctly when transforming from lat/lon to meter and vice-versa? Are
there any other methods besides the ones I mentioned?

Cheers
Daniel

On Thu, Oct 23, 2008 at 8:16 AM, Nikos Alexandris

<nikos.alexandris@felis.uni-freiburg.de> wrote:

On Thu, 2008-10-23 at 10:40 +0100, Corrado wrote:

Thanks Nikos, thanks Daniel!

Unfortunately, my map is a raster map, not a vector .... can I still
use
v.in.region and v.proj?

Thanks

[...]

Hi Corrado!

For sure!! You can use v.in.region anytime since it deals with the
current region and not with a specific map.

Regards, Nikos

Corrado wrote:

I have tried everything I could. I read the GRASS books chapters and the
documentation. It does not work.

The problem is that the files I have are in WGS84 unprojected. See here:

http://www.worldclim.org/format.htm

Whichever procedure I use, whther Daniel's or the book (for example: page 47),
the problem is always the same.

I create a location + mapset + map with the worldclim data (location01). I
load the original data, and I call the raster map "precipitation". I generate
the correct box with v.region and I call it "test". I create a new location +
mapset (utm + osgb), called location02. I run v.proj in location B pointing
at map "test"

v.proj input=test location=location01 mapset=PERMANENT

I have the same problem:

PROJ_INFO file not found for location location01

But of course the projection is not defined, the data are unprojected!

In this context, "lat/lon" is considered a projection. "unprojected"
or "X/Y" means that the coordinate system is entirely arbitrary (e.g.
pixel coordinates in an aerial photograph).

I am slowly getting desperate .... I have been working on it several days and
cannot get a sensible result. I have tried with different datum
transformation parameters for osgb36 and trying to project .... but the
results are very unreliable ....

Do you know what should I do?

Use g.proj or g.setproj to set the projection information for
location01, choosing lat/lon WGS84.

--
Glynn Clements <glynn@gclements.plus.com>

Corrado,
What v.in.region does is create a vector that surrounds your current
region. Then, when you run v.proj, the vector covering the region is
brought to your target location so you know where your projected
raster should be at...

The "manual labor" after that is setting the resolution correctly. You
could then use some rule of thumb like 3arcsec ~= 90 m or, try to set
the number of columns/rows in your target region similar to your
lat/lon region. But this latter option does not work so well if your
raster is to be tilted a lot...

On an "additional" question: How do people set the resolution
correctly when transforming from lat/lon to meter and vice-versa? Are
there any other methods besides the ones I mentioned?

Cheers
Daniel

On Thu, Oct 23, 2008 at 8:16 AM, Nikos Alexandris
<nikos.alexandris@felis.uni-freiburg.de> wrote:

On Thu, 2008-10-23 at 10:40 +0100, Corrado wrote:

Thanks Nikos, thanks Daniel!

Unfortunately, my map is a raster map, not a vector .... can I still
use
v.in.region and v.proj?

Thanks

[...]

Hi Corrado!

For sure!! You can use v.in.region anytime since it deals with the
current region and not with a specific map.

Regards, Nikos

Cheers,

--
Sebastian P. Luque
Department of Biology
Memorial University
St. John's, NL A1C 3X9
Canada
http://sebmags.homelinux.org/~sluque

Currently at:
Fisheries and Oceans Canada
501 University Crescent
Winnipeg, MB R3T 2N6
Canada
Tel +1 (204) 586-8170
Fax ... 984-2403

On Thu, 23 Oct 2008 08:38:47 -0200,
"Daniel Victoria" <daniel.victoria@gmail.com> wrote:

[...]

On an "additional" question: How do people set the resolution
correctly when transforming from lat/lon to meter and vice-versa? Are
there any other methods besides the ones I mentioned?

[Apologies for my previous botched post (got jumpy fingers today).]

I'm also looking for an optimized procedure to do this, so would
appreciate any followups to this question. Thanks.

--
Seb

Sebastian P. Luque wrote:

On an "additional" question: How do people set the resolution
correctly when transforming from lat/lon to meter and vice-versa?

There is no "correct" resolution when re-projecting rasters; whatever
you choose, the conversion will be "lossy" to some extent. It's just a
question of what's adequate for your purposes.

--
Glynn Clements <glynn@gclements.plus.com>

Alex Mandel wrote:

Even though the projection is not defined, the data still has some
information it is WGS84 which has a datum/coordinate system so you do
need to define it, the EPSG code is 4326.

You can't reproject something from nothing, it needs at least a
coordinate system, otherwise it doesn't know what point on your current
map should be moved and to where it should be moved.

When you configured your first location you should have specified the
EPSG as 4326.

right,

"They are in geodetic coordinate system (not projected, i.e., 'GEOGRAPHIC'
or 'LATLONG' system). The datum is WGS84."
   -- http://www.worldclim.org/format.htm

by "projected" they mean "not a grid system", they do not mean "without
specific georeferencing".

technically speaking lat/lon is not actually "projected" onto any plane.
it is just that it has come into common language to call any SRS a
projection, even when it is lat/lon.

so Alex is right, just use plain old lat/lon WGS84 epsg:4326 and be happy
it is not like Google Maps mercator-on-perfect-sphere.

Hamish

Dear friends,

this is what I have done:

1) created a location test12, coordinate system “latitude – longitude”, wgs84
datum, default datum transformation parameters (1), region 90N,60S
(=-60) ,180W (=-180),180E, resolution 30 seconds (=0.00833333 degrees) on
both the X and y axis.

2) imported the file: r.in.bin -s bytes=2 input=prec/prec_4.bil north=90N
south=60S west=180W east=180E output=dummy01 rows=18000 cols=43200
anull=-9999

3) selected the region of interest: g.region n=64N s=49N w=10W e=2E

4) calculated the width of the 30 seconds raster pixel, using the command
g.region -m:

<snippet>
n=64
s=49
w=-10
e=2
nsres=927.90359439
ewres=508.14568663
rows=1800
cols=1440
</snippet>

5) created the region box, using the command:

v.in.region output=region_box

6) created a new location, test13, using UTM, osgb36 datum, (1) datum
transformation parameters (5 m resolution), zone 30, region 1300000
(N),-200000 (S), -800000(W),800000(E), resolution 1 m.

7) loaded the region box:

v.proj input=region_box location=test12 mapset=PERMANENT
output=region_box_fromtest12

8) set the region to the region identified by the vector:

g.region vect=region_box_fromtest12

9) set the resolution to the resolution calculated previously:

g.region nsres=927.90359439 ewres=508.14568663

10) project the raster file:

r.proj input=dummy01 location=test12 mapset=PERMANENT
output=test_fromtest12_02 method=cubic

11) set the resolution I want:

g.region nsres=2000 ewres=2000

12) resample:
r.resample input=test_fromtest12_02
output=test_fromtest12_02_resample --overwrite

What do you think?

On Saturday 01 November 2008 04:33:55 Glynn Clements wrote:

Sebastian P. Luque wrote:
> On an "additional" question: How do people set the resolution
> correctly when transforming from lat/lon to meter and vice-versa?

There is no "correct" resolution when re-projecting rasters; whatever
you choose, the conversion will be "lossy" to some extent. It's just a
question of what's adequate for your purposes.

--
Corrado Topi

Global Climate Change & Biodiversity Indicators
Area 18,Department of Biology
University of York, York, YO10 5YW, UK
Phone: + 44 (0) 1904 328645, E-mail: ct529@york.ac.uk

Corrado wrote:

this is what I have done:

9) set the resolution to the resolution calculated previously:

g.region nsres=927.90359439 ewres=508.14568663

10) project the raster file:

r.proj input=dummy01 location=test12 mapset=PERMANENT
output=test_fromtest12_02 method=cubic

11) set the resolution I want:

g.region nsres=2000 ewres=2000

12) resample:
r.resample input=test_fromtest12_02
output=test_fromtest12_02_resample --overwrite

What do you think?

You would probably get better results if you replaced the above four
steps with just:

  g.region nsres=2000 ewres=2000
  r.proj ...

Resampling is invariably lossy, so fewer resampling steps are
preferable to more steps.

--
Glynn Clements <glynn@gclements.plus.com>

Dear GLynn,

does the r.proj resample the raster?

The documentation seems to suggest you have to use g.region -m, .... or may be
I did not understand it properly.

Here it is http://grass.osgeo.org/grass62/manuals/html62_user/r.proj.html, in
the section NOTES.

What do you think?

On Monday 03 November 2008 19:55:03 Glynn Clements wrote:

Corrado wrote:
> this is what I have done:
>
> 9) set the resolution to the resolution calculated previously:
>
> g.region nsres=927.90359439 ewres=508.14568663
>
> 10) project the raster file:
>
> r.proj input=dummy01 location=test12 mapset=PERMANENT
> output=test_fromtest12_02 method=cubic
>
> 11) set the resolution I want:
>
> g.region nsres=2000 ewres=2000
>
> 12) resample:
> r.resample input=test_fromtest12_02
> output=test_fromtest12_02_resample --overwrite
>
> What do you think?

You would probably get better results if you replaced the above four
steps with just:

  g.region nsres=2000 ewres=2000
  r.proj ...

Resampling is invariably lossy, so fewer resampling steps are
preferable to more steps.

--
Corrado Topi

Global Climate Change & Biodiversity Indicators
Area 18,Department of Biology
University of York, York, YO10 5YW, UK
Phone: + 44 (0) 1904 328645, E-mail: ct529@york.ac.uk

Hi List,

I’m trying to display a Quickbird imagery with nice real colors.

I tried with:

r.composite red=m1.red green=m1.green blue=m1.blue out=composite
i.fusion.brovey ms1=m1.green ms2=m1.nir ms3=m1.red pan=pan outputprefix=fusion -q

I tried also to set the RGB maps to grey und grey.eq but unfortunately until now the results are not very satisfying.

Any suggestion will be appreciated.

regards,
tommaso