4

I am attempting to use the point sampling plugin in QGIS to extract values from phenology layers for the Western US provided by USGS. The data comes a .hdr file with a prj file specifying a Lambert Azimuthal Equal-area projection.

PROJCS["Lambert Azimuthal Equal area",
GEOGCS["GCS_Sphere_ARC_INFO",DATUM["D_Sphere_ARC_INFO",SPHEROID["Sphere_ARC_INFO",6370997.0,0.0]],
PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
PROJECTION["Lambert_Azimuthal_Equal_Area"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",-100.0],
PARAMETER["Latitude_Of_Origin",45.0],
UNIT["Meter",1.0]]

Points for my sites were generated from a csv with lat long coords in decimal degree. I have tried placing the points into the same crs as the raster, I have tried modifying the false northing of the raster to align the points and the raster, I have tried reprojecting both to a third projection and modifying the false northing of the raster to make them align.

When I use a categorical raster with the same crs args as the point file, the extractions work as expected.

The result with the phenology raster extraction is empty cells in datatable.

Any ideas?

screen shot of aligned layers

(ignore the fact that the phenology data is likely highly correlated with the landcover data - boss wants it)

  • currently looking at USGS report on LAEA type projections (page 44). – Collin Feb 20 '17 at 18:22
3

Point sampling tool cannot work properly on data with different projections. The output will always be empty if the raster file and point shapefile have different projection. Since you tried to reproject the point data to have similar projection with the raster data.

One possible solution is to load the data on a new fresh QGIS project file, both the raster data and the projected point data, which have a similar projection, and try to use point sampling tool.

  • When I load the files in a fresh project, the result is all the points appear in a single cell of the raster. The coords for the shapefile appear are in meters but have the same values as the decimal degrees. When I change the CRS in R using the rgdal sptransform function, the points are north of the raster image. I think I am not understanding how the Central Meridian and latitude of origin parameters are impacting the outputs I see. – Collin Feb 20 '17 at 18:07
0

@ahmadhanb is correct that the layers must be aligned. I did not explicitly state that I understood that in my question and appreciate the clarity. To properly align the layers, I modified the crs string as follows:

original_crs<- +proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs

crs_cnvrt <- "+proj=laea +lat_0=45 +lon_0=-100 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs"

The lat_0 and lon_0 parameters were distorting the projection. By setting lat_0 to the latitude of origin and the lon_0 to the meridian of origin the projection worked.

In R, I used the rgdal::sptransform function to transform the points then exported the points as a csv. When uploading the points in QGIS I used the original CRS arguments from the raster layer when prompted.

#get sites
coords <- sites[,match(x = c("Latitude","Longitude"),names(sites))]

#Create spatial object and set CRS
names(coords) <- c("x", "y") #rename

coordinates(coords) <- ~ y + x

proj4string(coords) <- CRS("+proj=longlat +datum=WGS84 +no_defs")

#Project to laea 
crs_cnvrt <- "+proj=laea +lat_0=45 +lon_0=-100 +x_0=45 +y_0=0 +datum=WGS84 +units=m +no_defs"

sites_transformed <- spTransform(coords, CRS(crs_cnvrt)) #reprojecting to raster projection

#Create CSV to export to QGIS
sites_t <- as.data.frame(sites_transformed@coords)
sites_t$name <- sites$Site

write.csv(sites_t,file = "./SITES_forextract45.csv")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.