[GRASS-user] reprojection problem

Tyler_Smith2 · December 15, 2011, 5:11pm

Hi,

I’m a GIS novice, trying to reproject a vector layer. My objective is to get a world map
projected into Lambert Equal-Area, centered on the north pole. I chose this projection
as it shows the entire arctic well, but I don’t have to use Lambert Equal-Area
specifically.

This is what I’ve tried so far:

Created a new location with the location wizard, selecting
‘read projection and datum terms from georeferenced data file’

The target file was:
WorldCountries.shp

retrieved from http://www.sharegeo.ac.uk/handle/10672/124

This appears to be correct, as I now have my location world-ll, which contains in the
PERMANENT mapset the layer WorldCountries. The region info is:

g.region -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 83:37:38.7084N
south: 90S
west: 180W
east: 180E
nsres: 8:40:52.93542
ewres: 18
rows: 20
cols: 20
cells: 400

I then made a new location, using the wizard again:
select coordinate system from a list: laea
central parallel: 90
central meridian: 100
false northing: 0
false easting: 0
datum code: wgs84

I then set the extent to 90N, -90S, 180E, -180W, NS resolution 0.1, EW resolution 0.1

The resulting region info is:

g.proj -p
-PROJ_INFO-------------------------------------------------
name : Lambert Azimuthal Equal Area
proj : laea
datum : wgs84
ellps : wgs84
lat_0 : 90
lon_0 : 100
x_0 : 0
y_0 : 0
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : meter
units : meters
meters : 1

From here I used the gui to import the WorldCountries vector layer from world-ll,
generating the command:

v.proj input=WorldCountries location=world-ll mapset=PERMANENT output=WorldCountries

This command runs without error, but with a few suspicious warnings:

v.proj input=WorldCountries location=world-ll mapset=PERMANENT output=WorldCountries
Building topology for vector map …
Registering primitives…
55402 primitives registered
417249 vertices registered
Building areas…
Area of size = 0.0 ignored
Area of size = 0.0 ignored
27423 areas built
26809 isles built
Attaching islands…
Attaching centroids…
Number of nodes: 54788
Number of primitives: 55402
Number of points: 0
Number of lines: 0
Number of boundaries: 28000
Number of centroids: 27402
Number of areas: 27423
Number of isles: 26809
Number of incorrect boundaries: 1
Number of centroids outside area: 1
Number of areas without centroid: 22

At this point, the WorldCountries layer appears in the Map Layers tab, but not on the
display. Zooming to selected layer produces an error:
Unable to zoom to vector map WorldCountries@PERMANENT

What have I done wrong?

Thanks!

Tyler

Markus_GRASS · December 15, 2011, 8:05pm

On Thu, Dec 15, 2011 at 6:11 PM, Tyler Smith <tyler.smith@mail.mcgill.ca> wrote:

Hi,

I'm a GIS novice, trying to reproject a vector layer. My objective is to get
a world map
projected into Lambert Equal-Area, centered on the north pole. I chose this
projection
as it shows the entire arctic well, but I don't have to use Lambert
Equal-Area
specifically.

This is what I've tried so far:

Created a new location with the location wizard, selecting
'read projection and datum terms from georeferenced data file'

The target file was:
WorldCountries.shp

retrieved from http://www.sharegeo.ac.uk/handle/10672/124

This appears to be correct, as I now have my location world-ll, which
contains in the
PERMANENT mapset the layer WorldCountries. The region info is:

g.region
-p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 83:37:38.7084N
south: 90S
west: 180W
east: 180E
nsres: 8:40:52.93542
ewres: 18
rows: 20
cols: 20
cells: 400

I then made a new location, using the wizard again:
select coordinate system from a list: laea
central parallel: 90
central meridian: 100
false northing: 0
false easting: 0
datum code: wgs84

I then set the extent to 90N, -90S, 180E, -180W, NS resolution 0.1, EW
resolution 0.1

The resulting region info is:

g.proj
-p
-PROJ_INFO-------------------------------------------------
name : Lambert Azimuthal Equal Area
proj : laea
datum : wgs84
ellps : wgs84
lat_0 : 90
lon_0 : 100
x_0 : 0
y_0 : 0
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : meter
units : meters
meters : 1

From here I used the gui to import the WorldCountries vector layer from
world-ll,
generating the command:

v.proj input=WorldCountries location=world-ll mapset=PERMANENT
output=WorldCountries

This command runs without error, but with a few suspicious warnings:

v.proj input=WorldCountries location=world-ll mapset=PERMANENT
output=WorldCountries
Building topology for vector map <WorldCountries>...
Registering primitives...
55402 primitives registered
417249 vertices registered
Building areas...
Area of size = 0.0 ignored
Area of size = 0.0 ignored
27423 areas built
26809 isles built
Attaching islands...
Attaching centroids...
Number of nodes: 54788
Number of primitives: 55402
Number of points: 0
Number of lines: 0
Number of boundaries: 28000
Number of centroids: 27402
Number of areas: 27423
Number of isles: 26809
Number of incorrect boundaries: 1
Number of centroids outside area: 1
Number of areas without centroid: 22

At this point, the WorldCountries layer appears in the Map Layers tab, but
not on the
display. Zooming to selected layer produces an error:
Unable to zoom to vector map <WorldCountries@PERMANENT>

What have I done wrong?

In short, Lambert Azimuthal Equal Area projection centered on the
north pole can not represent the whole earth.

When I reproject the full vector, I get with v.info:

v.info -g map=WorldCountries@PERMANENT
north=inf
south=-12469749.9061447
east=inf
west=-12555936.84435579
top=0.000000
bottom=0.000000

Note the inf for north and easth. Something went wrong there, inf
means infinity.

Then I tried to import only the northern hemisphere because a Lambert
Azimuthal Equal Area projection is, AFAICT, not able to represent the
whole earth, particularly if this projection is centered on the north
pole.

g.region n=90 s=0 e=180 w=-180 -p

v.in.region output=northern_hemisphere
v.select ainput=WorldCountries@PERMANENT
binput=northern_hemisphere@PERMANENT output=WorldCountries_North

Don't be afraid of the warnings produced by v.select, they are harmless.

I closed GRASS and started it in the location defined by the Lambert
Azimuthal Equal Area projection centered on the north pole and
reprojected the country boundaries of the northern hemisphere with:

v.proj input=WorldCountries_North location=world_ll mapset=PERMANENT
output=WorldCountries

and everything (output of v.proj, output of v.info, display) is fine

Essentially, you must avoid Antactica which is including the south
pole, which in turn conflicts fatally with a Lambert Azimuthal Equal
Area projection centered on the north pole.

Markus M

Hamish · December 16, 2011, 1:40am

Markus Metz wrote:

In short, Lambert Azimuthal Equal Area projection centered
on the north pole can not represent the whole earth.

If you're focused on the N pole, Polar Stereographic projection may be
the way to go.

e.g. NASA's Blue Marble 2002 and Night on Earth reprojected from lat/long:
http://grass.osgeo.org/screenshots/images/day_on_earth_N.png
http://grass.osgeo.org/screenshots/images/night_on_earth_N.png

IIRC grid was done with with v.mkgrid in lat/long location then
reprojected, 'd.grid -g' isn't as robust.

(from http://grass.osgeo.org/screenshots/raster.php)

there are a few versions to choose from:
$ proj -lp | grep -i Stereo

2c,
Hamish