I am still having problems with getting the raster to scale and
match/overlay perfectly with my other datasets. To recap, I was able to
convert the netcdf to a raster format successfully, but the projections were
skewed. Robert suggested to include the proj string and hence define the
projection for the raster. Even after doing that, it fails to achieve the
desired result which I think is a consequence of the extents not being
specified explicitly/correctly. I am therefore including a header file entry
here, one which was created by a GIS expert in our research group who has
since left. I know it includes the extents but I am finding it a bit hard to
decipher. Here it is:
ENVI
description = {
File Imported into ENVI.}
samples = 507
lines = 356
bands = 1
header offset = 0
file type = ENVI Standard
data type = 5
interleave = bsq
sensor type = Unknown
byte order = 0
map info = {vulcan, 1.0000, 1.0000, -2736000.0348, 1952000.0061,
1.0000000000e+004, 1.0000000000e+004, WGS-84, units=Meters}
projection info = {4, 6378137.0, 6356752.3, 40.000000, -97.000000, 0.0, 0.0,
33.000000, 45.000000, WGS-84, vulcan, units=Meters}
coordinate system string =
{PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-97.0],PARAMETER["Standard_Parallel_1",33.0],PARAMETER["Standard_Parallel_2",45.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",40.0],UNIT["Meter",1.0]]}
wavelength units = Unknown
Does anyone know how I can go about setting the extents of the raster from
the above header file?
Thanks much,
Advait
On Sat, Nov 6, 2010 at 9:03 PM, Robert J. Hijmans <r.hijm...@gmail.com>
wrote:
> Advait,
>
> I overlooked the datum bit. You can try this:
>
> projection(r) <- "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-96
> +datum=NAD83"
>
> Robert
>
> On Sat, Nov 6, 2010 at 3:57 PM, Advait Godbole <advaitgodb...@gmail.com>
> wrote:
> > I was using Mike's "writeGDAL" function to output the file. USing
> > writeRaster as you have shown gives the following output:
> >> x = writeRaster(r, "monthlyHDR.10k.2002.annsum2.tif", overwrite = TRUE)
> >> x
> > class : RasterLayer
> > filename :
> >
> C:/Work/13Mar.research.bkup/MS_Thesis/work.data.analysis/spatial/climatezone/monthlyHDR.10k.2002.annsum2.tif
> > nrow : 356
> > ncol : 507
> > ncell : 180492
> > min value : 1199.185
> > max value : 321709.3
> > projection : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-96 +x_0=0
> > +y_0=0 +ellps=WGS84 +units=m +no_defs
> > xmin : -137.3312
> > xmax : -62.11259
> > ymin : 22.12693
> > ymax : 57.35954
> > xres : 0.1483602
> > yres : 0.09896801
> > I see that the "ellps" paramter is still WGS84. Moreover, I opened ArcMap
> > again and have attached a screenshot of the layer properties window.
> Here,
> > the datum comes out to be unknown which should be NAD1983. So I still
> have
> > to nudge it somehow to get there right?
> > Advait
> > On Sat, Nov 6, 2010 at 6:35 PM, Robert J. Hijmans <r.hijm...@gmail.com>
> > wrote:
> >>
> >> Are you sure? How do you know? Perhaps fooled by ArcMap (need to close
> >> and open it to get rid of its cache). Below I show that it works, in
> >> principle.
> >>
> >> 2sp stands for 2 standard parallels. There is another way to define
> >> the same projection with a single parallel.
> >>
> >> Robert
> >>
> >> r <- raster(system.file("external/test.grd", package="raster"))
> >> > r
> >> class : RasterLayer
> >> filename : C:/soft/R/R-2.12.0/library/raster/external/test.grd
> >> nrow : 115
> >> ncol : 80
> >> ncell : 9200
> >> min value : 128.434
> >> max value : 1805.78
> >> projection : +init=epsg:28992
> >> +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812
> >> +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889
> >> +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs
> >> extent : 178400, 181600, 329400, 334000 (xmin, xmax, ymin, ymax)
> >> resolution : 40, 40 (x, y)
> >>
> >> > projection(r) <- "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-96"
> >> > r
> >> class : RasterLayer
> >> filename : C:/soft/R/R-2.12.0/library/raster/external/test.grd
> >> nrow : 115
> >> ncol : 80
> >> ncell : 9200
> >> min value : 128.434
> >> max value : 1805.78
> >> projection : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-96
> >> +ellps=WGS84
> >> extent : 178400, 181600, 329400, 334000 (xmin, xmax, ymin, ymax)
> >> resolution : 40, 40 (x, y)
> >>
> >> > x = writeRaster(r, 'test.tif', overwrite=TRUE)
> >> > x
> >> class : RasterLayer
> >> filename : C:/Documents and Settings/rhijmans/My Documents/test.tif
> >> nrow : 115
> >> ncol : 80
> >> ncell : 9200
> >> min value : 128.434
> >> max value : 1805.78
> >> projection : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-96
> >> +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
> >> extent : 178400, 181600, 329400, 334000 (xmin, xmax, ymin, ymax)
> >> resolution : 40, 40 (x, y)
> >>
> >> > r = raster('test.tif')
> >> > projection(r)
> >> [1] "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0
> >> +ellps=WGS84 +units=m +no_defs"
> >> >
> >>
> >>
> >>
> >> On Sat, Nov 6, 2010 at 3:29 PM, Advait Godbole <advaitgodb...@gmail.com
> >
> >> wrote:
> >> > Dear Robert,
> >> > I just tried doing what you suggested but the geotiff is still written
> >> > with
> >> > WGS84. Probably the prj string needs to have an inclusion of the EPSG
> >> > code
> >> > somewhere - 9802 from the website you noted? By the way, what does
> "2SP"
> >> > in
> >> > the proj name denote?
> >> > Thanks,
> >> > Advait
>
On Sat, Nov 6, 2010 at 6:11 PM, Robert J. Hijmans <r.hijm...@gmail.com>wrote:
> Dear Advait,
>
> "netcdf does not really have a "projection" in GIS parlance". netcdf
> is flexible file format in which you can store all kinds of data in
> different ways. Climatologists and some others tend to follow the "CF"
> (NetCDF Climate and Forecast Metadata Convention) conventions.
> Unfortunately, they do not make use one of the standard projections
> description formats. I have not taken the time yet to parse the CF
> descriptions and into a proj.4 string. However, the projection
> metadata is stored in the raster object that you create from a nc (CF)
> file in the slot "prj"
>
> r <- raster('ncfile')
> r...@prj
>
> From that you can construct a proj.4, in your case probably:
>
> prj <- "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-96"
>
> (see
> http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp.html
> )
>
> and do:
>
> projection(r) <- prj
>
> and save the file.
>
> Robert
>
> On Sat, Nov 6, 2010 at 12:51 PM, Advait Godbole <advaitgodb...@gmail.com>
> wrote:
> > Ok, so I was able to convert and visualize everything but the projections
> > are skewed. Do I have to define a projection when reading it in via
> raster
> > or using writeGDAL? I have a custom co-ordinate system which is a
> modified
> > North American Lambert Conformal Conic (NAD 1983 datum) like so:
> > False_Easting: 0.00000000
> > False_Northing: 0.00000000
> > Central_Meridian: -96.00000000
> > Standard_Parallel_1: 33.00000000
> > Standard_Parallel_2: 45.00000000
> > Latitude_Of_Origin: 40.00000000
> > Linear Unit: Meter
> >
> > I would presume that writeGDAL takes care of the projection but my hunch
> > says otherwise since the netcdf does not really have a "projection" in
> GIS
> > parlance if I am not mistaken.
> >
> > Thanks!
> >
> >
> > On Sat, Nov 6, 2010 at 2:16 PM, Advait Godbole <advaitgodb...@gmail.com
> >wrote:
> >
> >> Thanks Mike and Sam. This worked perfectly - I am only now realizing how
> >> powerful R and these libraries are. Very elegant and smooth!
> >>
> >> cheers,
> >>
> >> adv
> >>
> >>
> >> On Sat, Nov 6, 2010 at 11:30 AM, Samuel Veloz <sdve...@ucdavis.edu>
> wrote:
> >>
> >>> You can skip the second step in Michael's code and just write to a tif
> >>> file using raster:
> >>>
> >>> r <- raster("ncfile.nc") #you may need to specify varname here (which
> >>> band you are using)
> >>> writeRaster(r,"r.tif")
> >>>
> >>> Sam
> >>>
> >>>
> >>> On Sat, Nov 6, 2010 at 7:19 AM, Advait Godbole <
> advaitgodb...@gmail.com>wrote:
> >>>
> >>>> Dear Members,
> >>>>
> >>>> I want to convert a 356 x 507 (lat,lon) 2D netcdf file into a raster
> so
> >>>> that I can read it in to image analysis/GIS programs. I initially
> adopted
> >>>> an
> >>>> approach whereby I first used NCL to convert it to a Fortran
> unformatted
> >>>> binary and then convert it to a formatted sequential format using
> >>>> Fortran.
> >>>> This was then read in to ENVI as a "flt" file. I would like to know if
> >>>> there
> >>>> is an elegant way in R to directly convert to a geospatial raster,
> more
> >>>> so
> >>>> because I recently switched to using a different cluster and am
> running
> >>>> into
> >>>> i/o problems with all those binary conversions.
> >>>>
> >>>> Regards,
> >>>>
> >>>> --
> >>>> Advait Godbole, BE Mech
> >>>> Doctoral Student, Global Carbon Cycle Lab
> >>>> Earth and Atmospheric Sciences
> >>>> Purdue University
> >>>> West Lafayette, IN 47906
> >>>> USA
> >>>> voice: +1 317 730 5235
> >>>>
> >>>> [[alternative HTML version deleted]]
> >>>>
> >>>> _______________________________________________
> >>>> R-sig-Geo mailing list
> >>>> R-sig-Geo@stat.math.ethz.ch
> >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>>>
> >>>
> >>>
> >>
> >>
> >> --
> >> advait godbole
> >>
> >
> >
> >
> > --
> > advait godbole
> >
> > [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo@stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
--
advait godbole
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
