[GRASS-user] reprojection problem

[Markus Metz]

On Thu, Dec 15, 2011 at 6:11 PM, Tyler Smith <tyler.smith at 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 at 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 at 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 at PERMANENT
binput=northern_hemisphere at 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