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

Paul Van Deusen Paul_VanDeusen at uml.edu
Fri Mar 16 10:18:35 EDT 2007


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 :-)
>>>>
>>>> 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-----
>>>
>>>       
>>
>>     
>
>   


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




More information about the grass-user mailing list