[GRASS-user] Alaska shapefile: Convert from -180,180 to 0-360

Fri Mar 16 08:56:32 EDT 2007

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 :-)
> >>
> >> Working in the maritime area around New Zealand, 180 is an issue I 
> >> deal with regularly :-(
> >>
> >> 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 at grass.itc.it
> >>> http://grass.itc.it/mailman/listinfo/grassuser
> >>
> >> _______________________________________________
> >> grassuser mailing list
> >> grassuser at 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 at nhh.no