[GRASS-user] Alaska shapefile crosses 180th parallel

GRASS GRASS user
1

I'm trying to read in the state boundary shapefile from http://nationalatlas.gov. However, the Aleutian Islands
cross the 180th parrallel, so longitude goes from +180 to -180 at that point. I end up with the end of the
Aleutian Islands getting disconnected. What is the solution? I used v.in.ogr and set up a location using the information that came
with the shapefile. See spatial information below..

Spatial_Reference_Information:
  Horizontal_Coordinate_System_Definition:
    Geographic:
      Latitude_Resolution: 0.000278
      Longitude_Resolution: 0.000278
      Geographic_Coordinate_Units: Decimal degrees
    Geodetic_Model:
      Horizontal_Datum_Name: North American Datum of 1983
      Ellipsoid_Name: GRS1980
      Semi-major_Axis: 6378137
      Denominator_of_Flattening_Ratio: 298.257222

2

Paul Van Deusen wrote:

I'm trying to read in the state boundary shapefile from http://nationalatlas.gov. However, the Aleutian Islands
cross the 180th parrallel, so longitude goes from +180 to -180 at that point. I end up with the end of the
Aleutian Islands getting disconnected. What is the solution? I used v.in.ogr and set up a location using the information that came
with the shapefile. See spatial information below.

Welcome to my world! (oceanic GIS around New Zealand :slight_smile:

Reproject to a non lat/long coordinate system, or convert the longitudes to a 0-360 space instead.

If you use PostGIS, shp2pgsql the shapefile, then run a query to translate all features with a component with long<0 by 360 on the X axis, then add these back to your original dataset, then pgsql2shp this to a new shapefile.

It works for me :slight_smile:

FYI, I'm submitting a request to EPSG to have 2 lat/long projections, so as well as EPSG:4326 which defaults to +-180, we have an EPSG code for the same projection with a 0-360 long space. I don't know how they'll feel about it, but we'll see...

Cheers

Brent Wood

Spatial_Reference_Information:
Horizontal_Coordinate_System_Definition:
   Geographic:
     Latitude_Resolution: 0.000278
     Longitude_Resolution: 0.000278
     Geographic_Coordinate_Units: Decimal degrees
   Geodetic_Model:
     Horizontal_Datum_Name: North American Datum of 1983
     Ellipsoid_Name: GRS1980
     Semi-major_Axis: 6378137
     Denominator_of_Flattening_Ratio: 298.257222

_______________________________________________
grassuser mailing list
grassuser@grass.itc.it
http://grass.itc.it/mailman/listinfo/grassuser

3

Paul Van Deusen wrote:

Welcome to my world :slight_smile:

Working in the maritime area around New Zealand, 180 is an issue I deal with regularly :frowning:

I suggest you add 360 to every longitude < 0. This moves the break to zero degrees instead of 180.

You can also use a variety of local projections to avoid the issue.

There is a parameter in the CVS version of proj4 to set the central meridian, defaults to zero, but a reprojection to WGS84 (EPSG:4326) with explicit parameters, & the central meridian of 18- works well....

Brent Wood

I'm trying to read in the state boundary shapefile from http://nationalatlas.gov. However, the Aleutian Islands
cross the 180th parrallel, so longitude goes from +180 to -180 at that point. I end up with the end of the
Aleutian Islands getting disconnected. What is the solution? I used v.in.ogr and set up a location using the information that came
with the shapefile. See spatial information below..

Spatial_Reference_Information:
Horizontal_Coordinate_System_Definition:
   Geographic:
     Latitude_Resolution: 0.000278
     Longitude_Resolution: 0.000278
     Geographic_Coordinate_Units: Decimal degrees
   Geodetic_Model:
     Horizontal_Datum_Name: North American Datum of 1983
     Ellipsoid_Name: GRS1980
     Semi-major_Axis: 6378137
     Denominator_of_Flattening_Ratio: 298.257222

_______________________________________________
grassuser mailing list
grassuser@grass.itc.it
http://grass.itc.it/mailman/listinfo/grassuser

4

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

If you are only interested in Alaska an perhaps part of Canada and Russia, you could always project the data to Alaska Albers:
- -PROJ_INFO-------------------------------------------------
name : Albers Equal Area
proj : aea
datum : nad83
a : 6378206.4
es : 0.006768657997291279
lat_1 : 55
lat_2 : 65
lat_0 : 50
lon_0 : -154
x_0 : 0
y_0 : 0
no_defs : defined
- -PROJ_UNITS------------------------------------------------
unit : meter
units : meters
meters : 1

On Mar 15, 2007, at 1:43 PM, Brent Wood wrote:

Paul Van Deusen wrote:

Welcome to my world :slight_smile:

Working in the maritime area around New Zealand, 180 is an issue I deal with regularly :frowning:

I suggest you add 360 to every longitude < 0. This moves the break to zero degrees instead of 180.

You can also use a variety of local projections to avoid the issue.

There is a parameter in the CVS version of proj4 to set the central meridian, defaults to zero, but a reprojection to WGS84 (EPSG:4326) with explicit parameters, & the central meridian of 18- works well....

Brent Wood

I'm trying to read in the state boundary shapefile from http://nationalatlas.gov. However, the Aleutian Islands
cross the 180th parrallel, so longitude goes from +180 to -180 at that point. I end up with the end of the
Aleutian Islands getting disconnected. What is the solution? I used v.in.ogr and set up a location using the information that came
with the shapefile. See spatial information below..

Spatial_Reference_Information:
Horizontal_Coordinate_System_Definition:
   Geographic:
     Latitude_Resolution: 0.000278
     Longitude_Resolution: 0.000278
     Geographic_Coordinate_Units: Decimal degrees
   Geodetic_Model:
     Horizontal_Datum_Name: North American Datum of 1983
     Ellipsoid_Name: GRS1980
     Semi-major_Axis: 6378137
     Denominator_of_Flattening_Ratio: 298.257222

_______________________________________________
grassuser mailing list
grassuser@grass.itc.it
http://grass.itc.it/mailman/listinfo/grassuser

_______________________________________________
grassuser mailing list
grassuser@grass.itc.it
http://grass.itc.it/mailman/listinfo/grassuser

- -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
Gary Sherman
Micro Resources: http://mrcc.com
   *Geospatial Hosting
   *Web Site Hosting
"We work virtually everywhere"
- -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.1 (Darwin)

iD8DBQFF+e+01zKuzV6goTgRAjjrAJ9Fp1YYqi03Ad44d2bKgMa39HTn2ACgpFPX
TuQxiQQJ+vIxZDwXwaNH/2s=
=qi5X
-----END PGP SIGNATURE-----

5

Thanks for the feedback. I ended up converting longitude to 0-360 as Brent suggested. Its fairly easy to
do with R. I included the little R script that I wrote in case anyone else might find it useful. In general, R
is good for manipulating shapefiles.

======R script follows============

##Use this script to pull a state out of the national shapefile that came
##from http://nationalatlas.gov. It also converts lon from
##-180,+180 to 0-360, and also demonstrates how to work with shapefiles

#boundingBox returns a matrix of points representing a bounding box of a submitted user polygon
boundingBox <- function(xy){
yMin <-min(xy[,2]); yMax <-max(xy[,2]);
xMin <-min(xy[,1]); xMax <-max(xy[,1]);
bb<-c(xMin,yMin, xMax, yMax);
return(bb);
}

library(maptools) #load the library
fName="statesp020"; #name of shapefile that needs sorting and trimming
fOut="02"; #output file name, AK fips code=02
ST="Alaska" #state name in dbf file

shp <- read.shape(paste(fName,".shp",sep="")) #read in the original shapefile
dbf <- read.dbf(paste(fName,".dbf",sep="")) #read in the original dbf file

cat("Wait while making a polyList \n")
shapes <- Map2poly(shp,region.id=NULL)

keep=dbf[,"STATE"]==ST

sub=subset(shapes,keep)

dbf2 = dbf[keep,]

#convert from +180,-180 to 0-360 (only needed for Alaska)
for(i in 1:length(sub)){
vec=sub[i][[1]][,1]<0;
sub[i][[1]][vec,1]= 360 + sub[i][[1]][vec,1]
attr(sub[i][[1]],"bbox") = boundingBox(sub[i][[1]])
centroid=attr(sub[i][[1]],"centroid")$x
if(centroid<0)attr(sub[i][[1]],"centroid")$x=360+centroid
}

write.polylistShape(sub,dbf2,fOut,factor2char = TRUE)

cat("\n its done \n")

============end of R script=================

Gary Sherman wrote:

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

If you are only interested in Alaska an perhaps part of Canada and Russia, you could always project the data to Alaska Albers:
- -PROJ_INFO-------------------------------------------------
name : Albers Equal Area
proj : aea
datum : nad83
a : 6378206.4
es : 0.006768657997291279
lat_1 : 55
lat_2 : 65
lat_0 : 50
lon_0 : -154
x_0 : 0
y_0 : 0
no_defs : defined
- -PROJ_UNITS------------------------------------------------
unit : meter
units : meters
meters : 1

On Mar 15, 2007, at 1:43 PM, Brent Wood wrote:

Paul Van Deusen wrote:

Welcome to my world :slight_smile:

Working in the maritime area around New Zealand, 180 is an issue I deal with regularly :frowning:

I suggest you add 360 to every longitude < 0. This moves the break to zero degrees instead of 180.

You can also use a variety of local projections to avoid the issue.

There is a parameter in the CVS version of proj4 to set the central meridian, defaults to zero, but a reprojection to WGS84 (EPSG:4326) with explicit parameters, & the central meridian of 18- works well....

Brent Wood

I'm trying to read in the state boundary shapefile from http://nationalatlas.gov. However, the Aleutian Islands
cross the 180th parrallel, so longitude goes from +180 to -180 at that point. I end up with the end of the
Aleutian Islands getting disconnected. What is the solution? I used v.in.ogr and set up a location using the information that came
with the shapefile. See spatial information below..

Spatial_Reference_Information:
Horizontal_Coordinate_System_Definition:
Geographic:
Latitude_Resolution: 0.000278
Longitude_Resolution: 0.000278
Geographic_Coordinate_Units: Decimal degrees
Geodetic_Model:
Horizontal_Datum_Name: North American Datum of 1983
Ellipsoid_Name: GRS1980
Semi-major_Axis: 6378137
Denominator_of_Flattening_Ratio: 298.257222

_______________________________________________
grassuser mailing list
grassuser@grass.itc.it
http://grass.itc.it/mailman/listinfo/grassuser

_______________________________________________
grassuser mailing list
grassuser@grass.itc.it
http://grass.itc.it/mailman/listinfo/grassuser

- -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
Gary Sherman
Micro Resources: http://mrcc.com
*Geospatial Hosting
*Web Site Hosting
"We work virtually everywhere"
- -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.1 (Darwin)

iD8DBQFF+e+01zKuzV6goTgRAjjrAJ9Fp1YYqi03Ad44d2bKgMa39HTn2ACgpFPX
TuQxiQQJ+vIxZDwXwaNH/2s=
=qi5X
-----END PGP SIGNATURE-----

--
Paul Van Deusen
NCASI
600 Suffolk St. Fifth Floor
Lowell, MA 01854
978-934-1948

6

On Fri, 16 Mar 2007, Paul Van Deusen wrote:

Thanks for the feedback. I ended up converting longitude to 0-360 as
Brent suggested. Its fairly easy to do with R. I included the little R
script that I wrote in case anyone else might find it useful. In
general, R is good for manipulating shapefiles.

Very nice! Was your result similar to what you would get from the
recenter() method in the sp package, where 360 degrees offset is the
default:

download:

http://www.census.gov/geo/cob/bdy/co/co00shp/co02_d00_shp.zip

and unzip - then start R in the working directory:

library(maptools)
co02 <- readShapePoly("co02_d00.shp", proj4string=CRS("+proj=longlat"))
plot(co02)
summary(co02)
co02a <- SpatialPolygonsDataFrame(recenter(co02), data=as(co02,
  "data.frame"))
plot(co02a)
summary(co02a)

These are using new-style classes rather than the ad-hoc old-style classes
returned by read.shape() directly, and can also be used with readOGR() in
the rgdal package to read many vector formats. But you are right that, in
R, you can get at both the attribute and position data directly if you
need to. Merging the multi-polygon counties would be a next step ...

Roger

======R script follows============

##Use this script to pull a state out of the national shapefile that came
##from http://nationalatlas.gov. It also converts lon from
##-180,+180 to 0-360, and also demonstrates how to work with shapefiles

#boundingBox returns a matrix of points representing a bounding box of a
submitted user polygon
boundingBox <- function(xy){
yMin <-min(xy[,2]); yMax <-max(xy[,2]);
xMin <-min(xy[,1]); xMax <-max(xy[,1]);
bb<-c(xMin,yMin, xMax, yMax);
return(bb);
}

library(maptools) #load the library
fName="statesp020"; #name of shapefile that needs sorting and trimming
fOut="02"; #output file name, AK fips code=02
ST="Alaska" #state name in dbf file

shp <- read.shape(paste(fName,".shp",sep="")) #read in the original
shapefile
dbf <- read.dbf(paste(fName,".dbf",sep="")) #read in the original dbf file

cat("Wait while making a polyList \n")
shapes <- Map2poly(shp,region.id=NULL)

keep=dbf[,"STATE"]==ST

sub=subset(shapes,keep)

dbf2 = dbf[keep,]

#convert from +180,-180 to 0-360 (only needed for Alaska)
for(i in 1:length(sub)){
vec=sub[i][[1]][,1]<0;
sub[i][[1]][vec,1]= 360 + sub[i][[1]][vec,1]
attr(sub[i][[1]],"bbox") = boundingBox(sub[i][[1]])
centroid=attr(sub[i][[1]],"centroid")$x
if(centroid<0)attr(sub[i][[1]],"centroid")$x=360+centroid
}

write.polylistShape(sub,dbf2,fOut,factor2char = TRUE)

cat("\n its done \n")

============end of R script=================

Gary Sherman wrote:
> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
> If you are only interested in Alaska an perhaps part of Canada and
> Russia, you could always project the data to Alaska Albers:
> - -PROJ_INFO-------------------------------------------------
> name : Albers Equal Area
> proj : aea
> datum : nad83
> a : 6378206.4
> es : 0.006768657997291279
> lat_1 : 55
> lat_2 : 65
> lat_0 : 50
> lon_0 : -154
> x_0 : 0
> y_0 : 0
> no_defs : defined
> - -PROJ_UNITS------------------------------------------------
> unit : meter
> units : meters
> meters : 1
>
>
> On Mar 15, 2007, at 1:43 PM, Brent Wood wrote:
>
>> Paul Van Deusen wrote:
>>
>> Welcome to my world :slight_smile:
>>
>> Working in the maritime area around New Zealand, 180 is an issue I
>> deal with regularly :frowning:
>>
>> I suggest you add 360 to every longitude < 0. This moves the break to
>> zero degrees instead of 180.
>>
>> You can also use a variety of local projections to avoid the issue.
>>
>> There is a parameter in the CVS version of proj4 to set the central
>> meridian, defaults to zero, but a reprojection to WGS84 (EPSG:4326)
>> with explicit parameters, & the central meridian of 18- works well....
>>
>> Brent Wood
>>
>>> I'm trying to read in the state boundary shapefile from
>>> http://nationalatlas.gov. However, the Aleutian Islands
>>> cross the 180th parrallel, so longitude goes from +180 to -180 at
>>> that point. I end up with the end of the
>>> Aleutian Islands getting disconnected. What is the solution? I used
>>> v.in.ogr and set up a location using the information that came
>>> with the shapefile. See spatial information below..
>>>
>>> Spatial_Reference_Information:
>>> Horizontal_Coordinate_System_Definition:
>>> Geographic:
>>> Latitude_Resolution: 0.000278
>>> Longitude_Resolution: 0.000278
>>> Geographic_Coordinate_Units: Decimal degrees
>>> Geodetic_Model:
>>> Horizontal_Datum_Name: North American Datum of 1983
>>> Ellipsoid_Name: GRS1980
>>> Semi-major_Axis: 6378137
>>> Denominator_of_Flattening_Ratio: 298.257222
>>>
>>> _______________________________________________
>>> grassuser mailing list
>>> grassuser@grass.itc.it
>>> http://grass.itc.it/mailman/listinfo/grassuser
>>
>> _______________________________________________
>> grassuser mailing list
>> grassuser@grass.itc.it
>> http://grass.itc.it/mailman/listinfo/grassuser
>
> - -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
> Gary Sherman
> Micro Resources: http://mrcc.com
> *Geospatial Hosting
> *Web Site Hosting
> "We work virtually everywhere"
> - -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
>
>
>
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.1 (Darwin)
>
> iD8DBQFF+e+01zKuzV6goTgRAjjrAJ9Fp1YYqi03Ad44d2bKgMa39HTn2ACgpFPX
> TuQxiQQJ+vIxZDwXwaNH/2s=
> =qi5X
> -----END PGP SIGNATURE-----
>

--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand@nhh.no

7

I didn't know about the recenter method, but I'm not surprised that R has it.
My results looked fine after I pulled the converted shapefile into grass.

I have used R many times to manipulate shapefiles. It looks like I need to read up on
the latest things in the maptools package.

Roger Bivand wrote:

On Fri, 16 Mar 2007, Paul Van Deusen wrote:

Thanks for the feedback. I ended up converting longitude to 0-360 as
Brent suggested. Its fairly easy to do with R. I included the little R
script that I wrote in case anyone else might find it useful. In
general, R is good for manipulating shapefiles.
    
Very nice! Was your result similar to what you would get from the recenter() method in the sp package, where 360 degrees offset is the default:

download:

http://www.census.gov/geo/cob/bdy/co/co00shp/co02_d00_shp.zip

and unzip - then start R in the working directory:

library(maptools)
co02 <- readShapePoly("co02_d00.shp", proj4string=CRS("+proj=longlat"))
plot(co02)
summary(co02)
co02a <- SpatialPolygonsDataFrame(recenter(co02), data=as(co02, "data.frame"))
plot(co02a)
summary(co02a)

These are using new-style classes rather than the ad-hoc old-style classes returned by read.shape() directly, and can also be used with readOGR() in the rgdal package to read many vector formats. But you are right that, in R, you can get at both the attribute and position data directly if you need to. Merging the multi-polygon counties would be a next step ...

Roger

======R script follows============

##Use this script to pull a state out of the national shapefile that came
##from http://nationalatlas.gov. It also converts lon from
##-180,+180 to 0-360, and also demonstrates how to work with shapefiles

#boundingBox returns a matrix of points representing a bounding box of a submitted user polygon
boundingBox <- function(xy){
yMin <-min(xy[,2]); yMax <-max(xy[,2]);
xMin <-min(xy[,1]); xMax <-max(xy[,1]);
bb<-c(xMin,yMin, xMax, yMax);
return(bb);
}

library(maptools) #load the library
fName="statesp020"; #name of shapefile that needs sorting and trimming
fOut="02"; #output file name, AK fips code=02
ST="Alaska" #state name in dbf file

shp <- read.shape(paste(fName,".shp",sep="")) #read in the original shapefile
dbf <- read.dbf(paste(fName,".dbf",sep="")) #read in the original dbf file

cat("Wait while making a polyList \n")
shapes <- Map2poly(shp,region.id=NULL)

keep=dbf[,"STATE"]==ST

sub=subset(shapes,keep)

dbf2 = dbf[keep,]

#convert from +180,-180 to 0-360 (only needed for Alaska)
for(i in 1:length(sub)){
vec=sub[i][[1]][,1]<0;
sub[i][[1]][vec,1]= 360 + sub[i][[1]][vec,1]
attr(sub[i][[1]],"bbox") = boundingBox(sub[i][[1]])
centroid=attr(sub[i][[1]],"centroid")$x
if(centroid<0)attr(sub[i][[1]],"centroid")$x=360+centroid
}

write.polylistShape(sub,dbf2,fOut,factor2char = TRUE)

cat("\n its done \n")

============end of R script=================

Gary Sherman wrote:
    

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

If you are only interested in Alaska an perhaps part of Canada and Russia, you could always project the data to Alaska Albers:
- -PROJ_INFO-------------------------------------------------
name : Albers Equal Area
proj : aea
datum : nad83
a : 6378206.4
es : 0.006768657997291279
lat_1 : 55
lat_2 : 65
lat_0 : 50
lon_0 : -154
x_0 : 0
y_0 : 0
no_defs : defined
- -PROJ_UNITS------------------------------------------------
unit : meter
units : meters
meters : 1

On Mar 15, 2007, at 1:43 PM, Brent Wood wrote:

Paul Van Deusen wrote:

Welcome to my world :slight_smile:

Working in the maritime area around New Zealand, 180 is an issue I deal with regularly :frowning:

I suggest you add 360 to every longitude < 0. This moves the break to zero degrees instead of 180.

You can also use a variety of local projections to avoid the issue.

There is a parameter in the CVS version of proj4 to set the central meridian, defaults to zero, but a reprojection to WGS84 (EPSG:4326) with explicit parameters, & the central meridian of 18- works well....

Brent Wood

I'm trying to read in the state boundary shapefile from http://nationalatlas.gov. However, the Aleutian Islands
cross the 180th parrallel, so longitude goes from +180 to -180 at that point. I end up with the end of the
Aleutian Islands getting disconnected. What is the solution? I used v.in.ogr and set up a location using the information that came
with the shapefile. See spatial information below..

Spatial_Reference_Information:
Horizontal_Coordinate_System_Definition:
Geographic:
Latitude_Resolution: 0.000278
Longitude_Resolution: 0.000278
Geographic_Coordinate_Units: Decimal degrees
Geodetic_Model:
Horizontal_Datum_Name: North American Datum of 1983
Ellipsoid_Name: GRS1980
Semi-major_Axis: 6378137
Denominator_of_Flattening_Ratio: 298.257222

_______________________________________________
grassuser mailing list
grassuser@grass.itc.it
http://grass.itc.it/mailman/listinfo/grassuser
          

_______________________________________________
grassuser mailing list
grassuser@grass.itc.it
http://grass.itc.it/mailman/listinfo/grassuser
        

- -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
Gary Sherman
Micro Resources: http://mrcc.com
*Geospatial Hosting
*Web Site Hosting
"We work virtually everywhere"
- -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.1 (Darwin)

iD8DBQFF+e+01zKuzV6goTgRAjjrAJ9Fp1YYqi03Ad44d2bKgMa39HTn2ACgpFPX
TuQxiQQJ+vIxZDwXwaNH/2s=
=qi5X
-----END PGP SIGNATURE-----

--
Paul Van Deusen
NCASI
600 Suffolk St. Fifth Floor
Lowell, MA 01854
978-934-1948