texts: Why does Merge ignore one specific shapefile?

texts:
I am working with ArcGIS 10.3.1.4959 (License type: Advanced), and am trying to merge 3 point shapefiles together using the Data Management Tools -> General -> Merge functionality from ArcToolbox.

Between the shapefiles the fields do not match. But what's important is that they share the same coordinate system and geometry, being all points.

When I apply the Merge tool, the final result is never consistent: sometimes it includes the last element of the 1st shapefile, and sometimes the 1st shapefile is not at all present. I always get all the elements in the last two shapefiles, however. All three shapefiles are within the same Geodatabase.

I have tried to use Data Management Tools -> General -> Append but I am forced to use the NO_TEST schema type, losing useful information.

Why does this happen, and what can I do to walk around this obstacle?


texts: I have to transform the list of polygons that is in vector form into a raster mask that should be geographically aligned with the NDVI image that I have. I was advised to use rasterio.features.rasterize. But I am not sure how to use the function.


texts: Getting .argmax() in a list appended to an image collection: Google Earth Engine

texts:
I have appended a list of elevation data to a filtered image collection in GEE as metadata. I now want to find the position of which elevation (e.g. 25 meters) has the highest pixel count using .argmax(). I can do this successfully for a singular image, but struggling to work out how this will work for an image collection as a whole (20 images +).
I can't share the entire script for work reasons. The bit I can show is:
var newVariable = MyCollection.map(function(imgcol){
  var histogramValues = ee.Image(imgcol).select('elevation').reduceRegion({
    reducer: ee.Reducer.autoHistogram(),
    geometry: geometry,
    scale: 2,
    maxPixels: 1e12,
  })
  return imgcol.set('histogramValues', histogramValues.values().get(0))
})
print('Histogram appended', autoHisto)**

This returns an image collection with appended histogram values (amount of pixels at a given elevation). I am now struggling to find out where the highest pixel count is.


texts: Editing and saving CSV file to use in plugin for QGIS?

texts:
For using in a plugin (PS Time Series Viewer), I want to edit a CSV file in the way plugin wants. Somehow I think there's a problem with the CSV file. I would like to tell you all my steps because I cannot diagnose the problem properly.

I edited with MS Excel and saved the CSV file. I called it into QGIS (Layer/Add Layer/Add delimited text layer). I saved as shapefile (export/save feature as) this point layer. At this step, when I clicking a point with PS Time Series tool, I get an error. The error message is "No time series values found for the selected point".

I attached my CSV file and the plugin's test data (It is given as shapefile, and when I try on this, the plugin works very well).

Original CSV file: https://www.dropbox.com/s/p5z529wwzzmu6cy/_MyOriginalCSV.csv?dl=0

Edited CSV file: https://www.dropbox.com/s/i67bex4kzxehmqw/_MyEditedCSV.csv?dl=0

Plugin's test data: https://www.dropbox.com/s/cs35nhwze77y86n/_PluginsTestShapeFile.rar?dl=0

Error message: https://www.dropbox.com/s/24c244qdl6565bj/Photo%2030.03.2019%2020%2007%2048.jpg?dl=0

Is there any problem with my CSV file or I missed another point?

QGIS version: 3.6
PS Time Series version: 0.3.1


texts: Proxy to work with OpenLayers and Geoserver Locally with IIS express

texts: I need to set a proxy to work with OpenLayers and Geoserver, which are both using a different port on my pc. I'm the site with openlayers is hosted IIS express and I'm not using Apache (so the proxy.cgi solution is not for me). I need it to perform a wmsgetfeatureinfo and, at this point, I really don't know what to put in OpenLayers.ProxyHost =. Can you help?


texts: Using gdaldem to match cumulative count colouring of QGIS?

texts:
I am trying to colour a single band grey raster (NDVI) using gdal to match the same colouring when I bring the raster into QGIS and colour it with a cumulative count cut (2%/98%). When I use the code:

  gdaldem color-relief raster_ndvi.tif colour_percent2test.txt -alpha  raster_ndvi_colour.tif

with my colour textfile:

98%  26   153  71
74%  161  209  106
50%  255  255  197
26%  254  177  101
2%   255   0    0
nv   0     0    0    0

I get a similar raster as the QGIS version, but I don't see as much red

This is the QGIS output:

And this is what I get from the gdal code:

I have also tried using a colour textfile:

100% 26   153  71
98%  26   153  71
74%  161  209  106
50%  255  255  197
26%  254  177  101
2%   255   0    0
0%   255   0    0
nv   0     0    0    0

and it ends up looking the exact same as the other attempt.

Any suggestions on how I can get my colouring to match QGIS without typing in the actual raster values in the colour ramp?

I want to be able to colour multiple rasters this way without needing to know the raster values to type them into the textfile each time (so 0.37 =  red, 0.504 = orange, etc.) as this would take more time and I may as well just use QGIS instead.


texts:
I have a shapefile layer whose CRS is WGS84. However, all points are off by some fraction of degree, roughly corresponding to a dozen of meters. Is there a way in QGIS to fix this offset?
I can more or less work out how much the offset is because I can relate the points to other layers whose alignment I trust more, but I don't know how to say to QGIS to apply that offset to the incorrect layer.
I thought about setting a custom CRS for the incorrect layer, but I couldn't find how to specify a WKT which basically says "like WGS84, but please add a certain fixed amount to all latitudes and all longitudes".


texts: Semi, third, quarter circle and other compound style symbols in QGIS?

texts:
I think I used quarter (and other) circles in QGIS in the past. But it might be as well that it was only in different software. Anyway now I am not able to find in symbols any option for compound (dispersed groups) style . Am I missing something?

If the only option is SVG symbols, do you know about some publicly available nice set of symbols to use for this?


texts: Convert polygons to lines ignoring redundant vertices in QGIS

texts:
Provided the following Polygons layer with 6 polygons (see image below), I would like to calculate each line's segment's angle (amongst other things), so I should convert those six polygons into lines, where each line should have only two vertices (start and end).

As far as I know (I am far from being an expert), I could do the following:

Convert polygons into lines with Polygons to lines algorithm, which results into 6 closed lines with the same attributes as the polygons in Polygon layer.
Split the lines into segments with GRASS' v.split and specify 2 as the maximum number of vertices in segments (see screenshot with the configuration I used) [EDIT: after  Spacedman comment, I have discovered that I could also use Explode lines algorithm]
[]4

That would work if it weren't that I happen to have a vertex with id 38 (see first screenshot) which is redundant, as it breaks one line that does not have any other intersections nor angle's changes into two segments.

Is there any way to deal with such redundant vertices? (Of course, the image is just a  sample of a much more complex layer with many more vertices like that and I cannot manually fix them).

I know I could use Simlify algorithm, but unfortunately, it would not solve my problems because currently only offers three ways of simplifiying:

Distance-based: I cannot use this method, as the offending vertex is in the middle of the segment and far beyond any reasonable threshold. If I used that distance as threshold I would dramatically change the geometry with much worse negative effects.
Snap to grid: I do not have any grid to snap vertices to
Area: I cannot use it, as I am working with lines.

Surprisingly there's no method to fix vertices within a line where its neighbouring segments have the same azimuth angle, which would solve my problem.

From what I read, v.clean provides a tool called bpol which may be what I am looking for. From its manual:

  bpol: break (topologically clean) polygons (imported from non
  topological format, like ShapeFile). Boundaries are broken on each
  point shared between 2 and more polygons where angles of segments are
  different

Unfortunately, I haven't succeeded in using it, and when I want to run the tool using Lines layer as input layer, I get unexpected results, as only 10 vertices are used as "slicing" vertices which result in 24 lines, many of them have several segments with different angles (see image below):

At this point I have two questions:

Is there any other way to achieve what I want to achieve?
Am I doing something wrong with v.clean (which, I have to admit that I have failed miserably on all my attempts to use v.clean inside QGIS 3.2 in Linux for other purposes)?


texts: Leaflet _onDrag throws TypeError: this._positions is undefined

texts:
I am creating a map using the MapQuest API with Leaflet.

Up until today the map worked fine in development and no errors were being shown in FF dev console.

This morning when I was giving a presentation of my progress of my work so far, I noticed errors started showing in the console - TypeError: this._positions is undefined.  The error points to leaflet.js:2716:4

    _onDrag : function() {
        if (this._map.options.inertia) {
            var t = this._lastTime = +new Date, e = this._lastPos = this._draggable._newPos;
            this._positions.push(e), this._times.push(t), t - this._times[0] > 200 && (this._positions.shift(), this._times.shift())
        }
        this._map.fire("move").fire("drag")
    },

This error is only thrown when I grab the map and drag it to see a location offscreen i.e. Monterey, CA is in the middle of the map, I grab the map at Monterey and drag up and left to put Soledad on the viewable portion of the map.  No errors are thrown when zooming in or out on the map using +/- or the mouse wheel, nor any when clicking on a marker to see what's in the popup.

Yesterday I did start adding some geoJSON functionality, but I couldn't get it to work.  After seeing these errors produced, I removed the geoJSON code thinking this was the culprit.  It was not.  The errors are still thrown.

The map still functions as needed, but I'm sure this will cause problems down the line.  Any ideas on what could be causing this?

My basic map code with some generated markers:


texts: Accuracy of ST_Contains on geometry column with WGS84 data

texts:
Apologies if this has been asked before; I found some similar questions (this mailing list post was particularly interesting: https://lists.osgeo.org/pipermail/postgis-users/2008-June/020353.html) but nothing quite satisfied my curiosity.

I have a table containing WGS84 data for a bunch of multipolygons (very large ones; multiple countries in most cases). I want to know which one contains an arbitrary point (where the point is also in WGS84 lat/long).

I am not sure if the geometry or the geography type is best suited for this. I implemented it first using geometry and it seems to work for some spot checks. I tried it with geography as well, and it was orders of magnitude slower.

1) Is geometry acceptable for modeling this, or will it have trouble?

2) If there will be some error using geometry, what are the error characteristics and magnitude to expect? (For example, will there be issues near boundaries of adjacent multipolygons, etc.)

EDIT: Since ST_Contains doesn't work with geography columns, I guess ST_Intersects would be the closest thing?


texts: Splitting overlapping polygons in a layer and creating new ones with merged attributes in R

texts:
I have a large shapefile with plenty of overlapping polygons. From the areas where they intersect, I would like to create new polygons (and thus split the original ones). I would like to also to merge the stored attributes - please note here I have field_ID numeric ID that would need to be stored either as a character split by comma or as a vector.
This should result in 6 polygons. A solution in R is preferable, but ArcMap would be good too. Thank you.

library(sp)
# Funciton creating a spatial polygon:
create_poly <- function(coordin,id_number){
  coordinates = matrix(c(coordin[1], coordin[3],
                         coordin[1], coordin[4],
                         coordin[2], coordin[4],
                         coordin[2], coordin[3],
                         coordin[1], coordin[3]),
                       ncol = 2, byrow = TRUE)
  P_1 = Polygon(coordinates)
  Ps_1 = SpatialPolygons(list(Polygons(list(P_1), ID = id_number)), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
  ID_1 <- sapply(slot(Ps_1, "polygons"), function(x) slot(x, "ID"))
  df <- data.frame(rep(0, length(ID_1)), row.names=ID_1)
  Ps1_df <- SpatialPolygonsDataFrame(Ps_1, df)
  Ps1_df
}

# Create shapefile with polygon boundaries of interest:
Pol_1 <- create_poly(c(80,100,20,40),1)
Pol_2 <- create_poly(c(90,120,30,50),1)
Pol_3 <- create_poly(c(130,150,30,50),1)
Pol_4 <- create_poly(c(145,170,40,60),1)

Pol_1@data$properties <- c("red")
Pol_2@data$properties <- c("blue")
Pol_3@data$properties <- c("green")
Pol_4@data$properties <- c("purple")

Pol_1@data$field_ID <- 1
Pol_2@data$field_ID <- 2
Pol_3@data$field_ID <- 3
Pol_4@data$field_ID <- 4

# Join the polygons into one layer:
Pol_large <- do.call(rbind, list(Pol_1, Pol_2, Pol_3, Pol_4))

plot(Pol_large, axes = TRUE)

I can intersect those polygons but in result I do get 8 - the intersected areas are double (as they have two different properties) and the original polygons are not cut.

library(raster)
inter1=raster::intersect(Pol_large, Pol_large)

I also tried changing the polygons into lines and polygonising again, but it does loose the attributes columns and would need additional processing:

library(rgeos)
pol_merge = gLineMerge(as(Pol_large,"SpatialLines"))
regions = gPolygonize(gNode(pol_merge))
#Show the 6 regions:
par(mar=c(0,0,0,0))
par(mfrow=c(2,3))
for(i in 1:length(regions)){plot(regions);plot(regions[i,],add=TRUE,col="red")}

gContains(Pol_large, regions, byid=TRUE)


texts:
I'm trying to import a ShapeFile inside my PostGIS database using the ogr2ogr, however I'm receiving an error:

  ogr2ogr: symbol lookup error: ogr2ogr: undefined symbol:
  _ZNK17GDALDriverManager14GetDriverCountEv

The command that I'm using is that:

  ogr2ogr -append -f "PostgreSQL" PG:"host=127.0.0.1
  dbname=DATABASE_NAME user=USER_NAME password=PASSWORD" points.shp
  -skipfailures

I'm looking for it in Google, but without success.

Would anyone have any idea what is happening?


texts: Using azimuthal orthographic projection in QGIS?

texts:
I'm trying to create a globe view in QGIS 2.6 with an Azimuthal Orthographic projection using NaturalEarth's 110m countries dataset.
I'm following Hamish Campbell's blog post instructions that have been replicated by other users (eg here).
Problem 1
Using Hamish's recommended projection parameters with my lats/longs: +proj=ortho +lat_0=36.5 +lon_0=-118.8 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs to define a custom projection gives this result:

Lots of artefacts - what am I doing wrong / is it a bug?.
Problem 2
But Hamish mentions the presence of artefacts (I wouldn;t have thought the sort of artefacts i'm encountering but...) and recommends using some python to clip the dataset to the circle of the globe from the projection view.
I won't reproduce the python Hamish publishes (you can view it in step 2 in his post linked above) but having saved the file as clipper.py in QGIS python path, imported into the QGIS python window (import clipper) and done clipper.doCall(iface,36.5,-118.8) I don't get a new layer generated in the map window.
Are there any obvious things I might be doing wrong?


texts: QGIS: Batch processing interface won't separate input files into distinct rows

texts:
I'm trying to analyze a folder-full of las files using functions in Whitebox Tools. I have no problem analyzing individual files; the issue has to deal with batch processing.
When I right click the tool in the processing toolbox and select >1 input file, the batch processing dialog concatenates the inputs into a single row, separating each with a ";" rather than breaking them out into separate rows. This causes the tool to fail, as it expects a single input.

This behavior is consistent for all the WBT tools, as well as graphical models I've created incorporating WBT tools. Note that the tools work fine on single input or if I manually add rows in the batch dialog using the '+' icon -- it's just specifying multiple inputs to the batch processing that is the problem. Also, it's not a behavior I see with other tools in different toolboxes (i.e., they populate distinct rows rather than concatenating inputs).

I could manually hit "+" and specify each file, but I have hundreds...

I'm using QGIS version 3.4.0-Madeira and the latest version of WBT from https://github.com/jblindsay/whitebox-tools#4-available-tools


texts:
I'm new to using the modeler function in QGIS.  I'm trying to follow the tutorial in the screenshot.  But when I follow the steps it ends up with my in/out algorithm not connected to the "counties layer/county name attribute" input with arrows.
How do I connect all these together?


texts: Adding files within file geodatabase (*.gdb) folder into ArcMap?

texts:
I'm trying to use ArcMap to make a basemap and I currently have an unzipped folder full of files including .atx, .gbdindexes, .gdbtable, .gdbtablx, and .spx. I'm well aware that these files cannot be 'added' into the ArcMap database as is. So I've made a .gdb file and I can find the .gdb in both the Catalog and 'Add Data' menus. However, I am not sure how to add the files mentioned above (atx, gdbindexes, etc) into the .gdb in order to add map data into the map.

Can someone help me out? I'm not familiar with these other file types and how the folders work together.


texts: Minimum mapping unit for raster on ArcGIS

texts:
So I have this huge national wetlands raster and I'm trying to set a minimum mapping unit of 5 acres. Meaning I want to delete clusters of pixels that are smaller than 5 acres. But the raster's attribute table only has one row - which means the entire raster is one massive clump of pixels. Is this possible? Or any suggestions on some other tool which might help me achieve that?

Cell size: 30x30
Number of bands: 1
Pixel type: signed integer
Pixel depth: 8 bit
Pyramids: level 9, resampling, nearest neighbor

Thanks!


texts: Why does the buffer take up whole screen?

texts:
I am trying to create buffers around the Hurricane Sandy path I created in QGIS (points are lat/long).

I have both points from each day she was tracted and the line I created from it with Points2One.

When I go to make buffers around either the line or the points, the buffer takes up the whole screen. Its all just one color.

I have made buffers before but from data I did not create.

Right now all of my projections are set to WGS 84 EPSG 4326.


texts: QGIS 3.6 QUICKOSM buildings query CRS with metre untis to create buffer

texts:
Using QGIS 3.6.2

I'm running a QUICKOSM query to find residential buildings in Scotland. I want to then put a buffer around them. The query is fine, if a little slow. The only CRS I can get to accurately project them (i.e. not of the african coast or the coast of Cornwall) are OSGB 1936:4277 and WGS 84:4326. Both these have units in degrees and therefore I can't apply a buffer. Other OSGB (27700 and 7405) don't project correctly nor does Pseudo-Mercator (the Project CRS). OTFP is on by default in 3.6.

Does anyone know a CRS for Scotland/UK that had m units?


texts: Labeling rule depending on whether label fits in polygon

texts:
I'm labeling some country polygons with their names. The data have one attribute with the full country name and another one with a short version (e.g. Sweden/SE). I want to label the countries with the long version if it fits completely within the polygon and with the short version otherwise. Is it possible to solve this with rule-based labeling in QGIS?

I know about the draw option "Only draw labels which fit completely within feature". I'm trying to pour this option into a rule if that makes sense.

QGIS version: 3.10.1-A Coruña


texts: Marker shows as a dot instead of the icon image in openlayers

texts:
i am trying to add icons on my map using features but instead of showing the icosn i get a dot on the map in that position.

let map = new Map({
target: "map",
layers:[
new TileLayer({
source: OSM(),
})
],
view: new View({
center: [0,0],
zoom:0
})
});

let a = new Feature({
    geometry: new Point([100,100]),
    name: 'a',
    imgSrc: "assets/arrow.png"
  })
let vectorSource = new sourceLayer({
  features: [a],
  wrapX: false,
  });
let markerVectorLayer = new vectorLayer({
    source: vectorSource,
 });
map.addLayer(markerVectorLayer);


texts:
In the query below I am exporting only the endpoint of linestrings in my table routes to a json file.

My objective, however, is to export to a csv, with four separate columns - starttime, endtime, lat and lon

ogr2ogr -f GeoJSON endpoints.json  "PG:host=localhost dbname=dbname user=user password=password"  -sql "SELECT ST_AsText(ST_EndPoint(path::geometry)) from routes"

How could I change the query in psql below to include latitude and longitude of the endpoints as separate columns, using ST_AsText(ST_EndPoint(path::geometry)) ?

   \COPY (
        SELECT starttime, endtime, lat, lon from routes
    ) TO 'endpoints.csv' WITH CSV HEADER DELIMITER ',';


texts: Exporting mxd to image with transparent background from ArcGIS for Desktop?

texts:
I cannot get ArcMap 10.0 or 10.1 to "Export Map" to a png or gif with a transparent background -- no matter what I set the background and transparent colors to.  Thus, I created the following script, which does the trick, but I'd rather be able to get the default tool in ArcMap to do it.  Any better ideas?

# Purpose: This script exports an MXD view to a png file with transparent background,
# which is problematic to do manually due to a Windows-related bug in the ArcMap UI.
# Warning: It will make any white features (e.g. the default background color) transparent.
#
# INSTRUCTIONS:
# When in ArcMap, viewing a saved mxd...
# 1. From the Geoprocessing menu, select 'Python'.
# 2. In the new Python window, right-click somewhere after >>>, and select 'Load...'
# 3. Select this script.
# 4. When this script appears, press Enter twice to run it.
# 5. A new [mxd-name].png with transparent background will appear in the same directory as the mxd.
# 6. When finished/satisfied, close the Python window.
# http://forums.arcgis.com/threads/65463-Zoom-to-tile-and-Export-to-JPEG-with-World-File
import arcpy
png_resolution = 300    # Can adjust DPI
# Set mxd and png name
mxd = arcpy.mapping.MapDocument("CURRENT")
pngname = mxd.filePath.replace(".mxd",".png")
df = arcpy.mapping.ListDataFrames(mxd)[0]
# Calculate ideal image size from the mxd view and desired DPI
png_width = int((df.extent.XMax - df.extent.XMin) * png_res * 12 / df.scale)
png_height = int((df.extent.YMax - df.extent.YMin) * png_res * 12 / df.scale)
#Export PNG
print "Exporting: {0}\n{1} x {2} pixels, {3}-dpi".format(pngname,png_width,png_height,png_res)
arcpy.mapping.ExportToPNG(mxd, pngname, df, df_export_width = png_width, df_export_height = png_height, resolution = png_res, transparent_color = "255, 255, 255")
#Clean up
del df, pngname, png_res, png_width, png_height, mxd


texts: Seeking free Web Processing Service (WPS) client?

texts:
A number of products seem to be implementing the server side of the OGC Web Processing Service (WPS 1.0.0) and these are listed at the OGC's Implementations by Specification page.

Is there a free WPS client software that can be installed on Windows and be ready to connect to a WPS 1.0.0 server from another vendor/team?


texts: Use harmonized world soil database (HWSD) in QGIS ie load .blw, .hdr and .stx files coming with .bil raster

texts:
I am trying to use the data of the HWSD in QGIS.
See: http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML/

How do you get the same result in QGIS as in the viewer ?

More precisely :
It is easy to open the .bil file in QGIS and open a table extracted from the database (HWSD_Q).
However, I can't figure out how to color each pixel to get the map you can see in the viewer.

One solution is to convert the raster to polygon (see extract data from HWSD in qgis) but this process is really too long and does not provide the expected result.

I guess the map is displayed correctly if the .blw, .hdr and .stx files are loaded in QGIS but is seems that these formats are not supported by QGIS. Futhermore, the "official" legend is very important to me as it shows the links between soils. Open the .bil raster and use Properties -> Style -> Singleband pseudocolor (or others) and classify is not a satisfying method because the process to manualy define the color ramp is too long.

I am using the latest LTR version of QGIS.


texts: How to generate maps like some at Wikipedia?

texts:
I've absolutely no idea if this is a proper forum to ask this or not. But I know that I can at least get some information here.

I was reading few wiki articles on Indian rivers and I'm very much impressed by the maps on those pages. Example:

Maps in the school geography textbook that I'm reading are dull & boring. Example:

These maps got me excited. I want to generate few maps like these. I specifically want to add details of the second map in to first map, for all the rivers shown in the first map. I also want to add names of reservoirs & dams with icons etc.

The guys who made these maps have also made lot of such maps for wiki:

http://commons.wikimedia.org/wiki/User:NordNordWest/Gallery/Maps
http://en.wikipedia.org/wiki/User:Shannon1/Maps

How can I generate such maps?

I'm an absolute newbie. I just discovered this.

Don't we have an advanced version of this tool? mapsofindia.com/editable-maps/customise-your-map.html where we can just select what we want to include in the map by using checkboxes


texts:
I'm having a similar problem to this post: Using Qgis server with Ubuntu 12.04

I am using Mint 15 desktop and QGIS 2.0.1.

I read the above post and decided to go back a step or two and followed this tutorial: http://docs.qgis.org/2.0/html/en/docs/user_manual/working_with_ogc/ogc_server_support.html

I only got to the testing part. I am able to connect to localhost/cgi-bin/qgis_mapserv.fcgi using the 'Add WMS/WMST Layers' tool and I can see 'wfslayer' in the list. The 'Add' button is greyed out and the dialogue says Coordinate Reference System (0 available).

When I create my own QGIS project the 'Add' button is also greyed out unless I add the shapefile's CRS to Project Properties>'CRS Exclusions' list. If I add the CRS then I can add the layer to the TOC but it does not render.

I started to check if my web server (apache2) configuration was correct as per: http://karlinapp.ethz.ch/qgis_wms/configuration/index.html but am having the following issues:

I have two fcgid.conf so do I edit both?
I do not know the 'syntax' of 'DefaultInitEnv LD_LIBRARY_PATH' for apache2.conf or fcgid.conf
I used /usr/lib/qgis for QGIS libraries (is this correct?) but cannot determine the path to 'qt libraries'.

Do I need to configure these 3 conf files?

Are they the likely cause of not being able to add the test WMS layer?


texts: Google Earth plugin deprecated. Which alternatives?

texts:
I've got a costumer of mine asking for an integration of the Google Earth Plugin to his website.

From Google earth Plugin site I read this:

  Google Earth Plugin API has been deprecated as of December 12th, 2014.
  The API will continue to work on supported browsers until December
  12th, 2015, and will shut down on that date.

What alternative exist?

I know that there are Google Maps API, Google Maps Engine API, Google Street View API, but I'm confused about them and which one can "substitute" Google Earth Plugin.

The features that I need are these, in order of importance:

load a KMZ file over a satellite view
show balloons with HTML data that is present in KMZ files entities. The html will contain links to street view, and I would like to have street view enabled inside balloons. I managed to do this in Google Earth.
satellite images historical imagery selector
measure tool
get address of a clicked point (geocooding)
3D buildings


texts: How to copy attributes after splitting polylines?

texts:
I have two almost identical shapefiles (polylines) representing a road network. The only difference is that the original shapefile has lines that have different lengths, while the other shapefile has a set length of 1000 meter.

While creating the equal sections of 1000 meters, the attribute table of the original shapefile with different lengths was lost. This also means that the new shapefile contains more Objects than the original shapefile.

Is there a way to copy the attributes from the original shapefile to the new shapefile? The attributes from the original need to be copied to all the objects that are at the same location in the new shapefile.

The first picture contains the original layer with classification as Major Road (please note, between two intersections there is only 1 object, unlike in the new layer). The second picture contains the new layer without classification and the third one contains the picture when using spatial join in ArcGis.


texts: how to divide the pyqgis code into many python modules

texts:
I created python script with one function that connect to mysql db, retrieves non spatial data, creates features and saves them into new layer.
Below is the pseuodo code:

###mysql_select.py
def create_features (mysql_connection)
    layer = QgsVectorLayer(xxx)
    provider = layer.dataProvider()
    QgsMapLayerRegistry.instance().addMapLayer(layer)
    fields = provider.fields()
    features = []

    mysql_connection.cursor.execute("select * from tables")

    for res1,res2,re3 in mysql_connection.cursor.fetchall():
    #store results to a layer attributes

    provider.addFeatures(features)
    _writer = QgsVectorFileWriter.writeAsVectorFormat(layer, path_to_shape_file)

###my_main_script.py
#import mysql_select
myConnection = MySQLdb.connect(db_credentials)
create_features(myConnection)
myConnection.close()

It works fine. But I want to divide all stuff to modules located in different .py files to get more comfortable development.
I modified the code above to mysql_select.py and my_main_script.py according to the marks ### above.
But I get the error "NameError: global name 'QgsVectorLayer' is not defined".

What is the problem? How to send layer variable to the external module?
What is the best method to divide the code into different modules in this case?
I have to make four sql selects and create 4 new layers, so to put all that stuff in one script seems not versatile to develop.

I`ve tried to put from qgis.core import QgsVectorLayer in the mysql_select.py and in the my_main_script.py but unsuccessful.


texts: Raster conversion from netCDF ORCA025 to equirectangular projection

texts:
I need to use raster data from Copernicus Marine Service, that comes as netCDF format, and with ORCA025 tri-polar grid (1440 x 1021 grid points, the three poles are located over Antarctic, Central Asia and North Canada).

But what I need it's just simple rasters using regular latitudes and longitudes (I think that's called equirectangular projection, WSG84...) I don't have a GIS software (and even, not enough knowledge to use them), just R.

So I guess that the conversion could be done with R and maybe gdal package, but I have no idea of how to do it.


texts: Why does the postGIS file lose geometry?

texts: I am loading shapefile data from a GPS into PostGIS via Leica ZenoOffice. This apparently imports fine with it's geometry in the attribute table. However, when I add a postGIS layer in QGIS the file uploads but without any geometry (all the other attributes have imported correctly). Is there something specific that this could be such as something simple like a box I have forgotten to tick or is this likely to be a problem with the GPS data itself?


texts: How to make points with label info in ArcGIS API for Silverlight?

texts:
I'm a new user of the Arcgis for Silverlight SDK.  Using c#, I would like to add points with labels showing the id for point.

How can I do this?


texts: Getting "iteration not started" when using UpdateCursor / deleteRow()?

texts:
I'm attempting to delete some features from a feature class using an arcpy.da.UpdateCursor. The line uc.deleteRow() throws an exception with the message iteration not started.

The where_clause is "item_no in (123,456,789)"; return records with a value in this field found in this comma-delimited list. Here is the where clause being passed in ITEM_NO in (1,8,9).

To remove the where_clause as a possible source of this problem, I hard-coded a single value into the where; "item_no = 123". This also threw the error.

If I change this to a searchCursor using the same parameters, three records are returned, as expected.

with arcpy.da.UpdateCursor(fc, "OBJECTID", where_clause=strWhere) as uc:
    for row in XYZ:
        try:
            uc.deleteRow()
        except Exception, e:
            print e

What am I doing wrong and how do I get around this?


texts: PyQGIS - get ellipsoid acronym of coordinate system of layer

texts:
I  want to calculate the area of polygons in a layer with PyQGIS.
I am using the following code:
d= QgsDistanceArea()

for feature in layer.getFeatures():
        area = d.measureArea(feature.geometry())
        print("Area : {}".format(area))

The Problem is, the coordinate sytem I use is ETRS89 / UTM zone 32 N (EPSG:25832). Therefore the returned values are not correct. As far as I know, to get correct values I have to set the ellipsoid (GRS 1980) like this:
d.setEllipsoid('ACRONYM_OF_ELLIPSOID')

What is the best way to get the acronym of the used ellipsoid?
Is it possible to read the acronym of the ellipsoid from the coordinate system? And if yes, how to do it?


texts: How to keep web map from displaying content from 2 zoom levels at the same time?

texts:
So I'm making a web app with a map that combines several layers (and dynamically toggles them off & on). I've made these layers in TileMill, exported them to .mbtiles files, and I'm serving them with TileStream. The problem I'm having is that these layers can be sparse -- there places that don't have any content. This leads to the map displaying content from two zoom levels at the same time.

The cause of this is that two well-meaning features that are interacting in an undesirable way:

TileMill doesn't generate tiles in areas where the layer has no content
When a tile is missing at one zoom level, Polymaps scales in a tile from another level. (So does Wax.)

Here's a picture of the problem. In this map, the street names and the (B) marker are on their own layer. Note the duplication for the (B) marker and the label of West 4th Ave. (in the lower left).

This diagram illustrates what's happening. The layer is blank in the upper left (since there's nothing to label in the water). Since TileMill didn't generate a tile there at zoom 15, Polymaps brings in the tile from zoom 14. If the tiles had opaque backgrounds, this would be OK. But since the tiles are transparent, the zoom 14 and 15 tiles overlap.

What can I do to avoid or work around this problem? I've come up with a few possibilities, all of which have drawbacks. Are there other options I'm overlooking? I'd appreciate any feedback or suggestions.

Switching toolkits: I'm open to moving to Leaflet or OpenLayers. But before I put in the effort, I'd appreciate some assurance that I won't have the same problem with those toolkits.

Using TileMill hosting: This will composite the tiles for me, but cost is a barrier. (I'd exceed the upload limits on the standard plan.) There's also the logistical complexity of making different composites for all the different combinations of layers I have.

Hacking up the .mbtiles files: The .mbtiles file is just a database, and I could add empty PNG files in for the missing tiles. This seems extreme, though, and I keep thinking there ought to be a simpler solution.


texts: Editing existing coordinate information with ArcPy?

texts:
I would like to find a way to edit shapefile (polyline) coordinate data using python. Basically, I want to edit all the z-coordinates in my shapefile so that they are 1 metre lower than they previously were.

I noticed that coordinate data can be manipulated en masse in the Esri Python tutorial video (https://www.youtube.com/watch?v=ge_U_w0PcFM, at 18:17 for anyone thats interested). They use a function within the field calculator to move all the x coordinates +5000m. I have reproduced the script below:

def shiftXCoord(Shape)
  point = shape.getpart(0)
  point.X  point.X + 5000
  return point
shiftXCoord(!Shape!)

I am experiencing two main problems using this script as a solution to my problem. (1) The commentators do not go into any great detail about what particular parts of the script mean/are used for. I know the above is a function but I do not know what lines 2 and 3 represent. (2) Also, without the field calculator I am unsure as to how to link this to the shapefile I am working on.

I am using ArcGIS Desktop 10.0.


texts: how to capture the clicked point before moving a vertex

texts:
I need the location of clicked point before moving a vertex.

IVertexMoved event will provide me the coordinates of the point after vertex is moved but I need before moving.

Please help.

Regards,
Sheethal


texts: Splitting line by other lines using ArcMap?

texts:
I have two shapefiles. The first: a line. The second: several lines overlapping the first line.

For some reason, I can't split the first line manually so now I want to split it where it's being overlapped by the other lines. Hopefully, this is possible without having to click on each overlapping feature in editor mode (I can't click on anything at all in editor mode).


texts:
How to show labels for geojson points in a Leaflet map?

There is Leaflet.label that is now deprecated in favour of L.Tooltip, but that only shows text on hover. I want to show the text labels directly on the map without needing user interaction.

Sample input - https://gist.github.com/maphew/e168430e999fc738eeb3448feda121cd

Desired result (green points with text labels, the other graphic elements are just for context):

Update: I want to create text that blends in with the map like the image below, not a popup tooltip.


texts: Creating lines connecting each possible pair of points in PostGIS

texts:
I am working with a group of points I would like to create lines connecting each possible pair of points.

Example: If I have points A, B and C I would like to have lines connecting A-B, A-C and B-C.

How can I do this in PostGIS?


texts: Using nearest neighbour with conditional field using PyQGIS?

texts:
I'm using QGIS 3.6

I have a point layer with the 84000 animal positions on a fence on varying dates.

I would like to obtain the nearest other animal to each animal point, on the condition the Date is the same. I have been able to add the nearest animal ID and its distance, but only by one date selection at a time.

Can I select all to do it at the same time, but on the condition the date matches?

I tried using a python code base on a different question but for QGIS 2.18 and I can't adapt it for more recent versions... maybe someone else can, link below for reference.

Finding closest segment to point which contains same attribute using NNJoin plugin of QGIS?


texts: How do I georeference automatically in QGIS 3.6

texts: I use QGIS 3.6 on Windows. I use HiRISE images and I want to co-register them to other background images. I use the georeferencer plugin and I do it manually but still, there is a pixel misalignment between these two images. Even though I do it many times it would be great if there is an automatic way to avoid this procedure.


texts: How to speed up the plotting of polygons in R?

texts:
I want to plot the country borders of North America over a raster image depicting some variable and then overlay contours on top of the plot using R. I have been successful in doing this using base graphics and lattice, but it seems that the plotting process is much too slow! I have not done this in ggplot2 yet, but I doubt that it will fare better in terms of speed.

I have the data in a netcdf file created from a grib file. For now, I downloaded the country borders for Canada, USA and Mexico, which were available in RData files from GADM which read into R as SpatialPolygonsDataFrame objects.

Here is some code:

# Load packages
library(raster)
#library(ncdf) # If you cannot install ncdf4
library(ncdf4)

# Read in the file, get the 13th layer
# fn <- 'path_to_file'
r <- raster(fn, band=13)

# Set the projection and extent
p4 <- "+proj=lcc +lat_1=50.0 +lat_2=50.0 +units=km +x_0=32.46341 +y_0=32.46341 +lon_0=-107 +lat_0=1.0"
projection(r) <- CRS(p4)
extent(r) <- c(-5648.71, 5680.72, 1481.40, 10430.62)

# Get the country borders
# This will download the RData files to your working directory
can<-getData('GADM', country="CAN", level=1)
usa<-getData('GADM', country="USA", level=1)
mex<-getData('GADM', country="MEX", level=1)

# Project to model grid
can_p <- spTransform(can, CRS(p4))
usa_p <- spTransform(usa, CRS(p4))
mex_p <- spTransform(mex, CRS(p4))

### USING BASE GRAPHICS
par(mar=c(0,0,0,0))
# Plot the raster
bins <- 100
plot(r, axes=FALSE, box=FALSE, legend=FALSE,
     col=rev( rainbow(bins,start=0,end=1) ),
     breaks=seq(4500,6000,length.out=bins))
plot(r, legend.only=TRUE, col=rev( rainbow(bins,start=0,end=1)),
     legend.width=0.5, legend.shrink=0.75,
     breaks=seq(4500,6000,length.out=bins),
     axis.args=list(at=seq(4500,6000,length.out=11),
                labels=seq(4500,6000,length.out=11),
                cex.axis=0.5),
     legend.args=list(text='Height (m)', side=4, font=2,
                      line=2, cex=0.8))
# Plot the borders
# These are so slow!!
plot(can_p, add=TRUE, border='white', lwd=2)
plot(usa_p, add=TRUE, border='white', lwd=2)
plot(mex_p, add=TRUE, border='white', lwd=2)
# Add the contours
contour(r, add=TRUE, nlevel=5)

### USING LATTICE
library(rasterVis)

# Some settings for our themes
myTheme <- RdBuTheme()
myTheme$axis.line$col<-"transparent"
myTheme$add.line$alpha <- 1
myTheme2 <- myTheme
myTheme2$regions$col <- 'transparent'
myTheme2$add.text$cex <- 0.7
myTheme2$add.line$lwd <- 1
myTheme2$add.line$alpha <- 0.8

# Get JUST the contour lines
contours <- contourplot(r, margin=FALSE, scales=list(draw=FALSE),
                        par.settings=myTheme2, pretty=TRUE, key=NULL, cuts=5,
                        labels=TRUE)

# Plot the colour
levels <- levelplot(r, contour=FALSE, margin=FALSE, scales=list(draw=FALSE),
                    par.settings = myTheme, cuts=100)

# Plot!
levels +
  layer(sp.polygons(can_p, col='green', lwd=2)) +
  layer(sp.polygons(usa_p, col='green', lwd=2)) +
  layer(sp.polygons(mex_p, col='green', lwd=2)) +
  contours

Is there a way to speed up the plotting of the polygons? On the system that I am working on, the plotting takes several minutes. I want to eventually make a function that will easily generate a number of these plots for inspection, and I reckon I will be plotting many of these maps, so I want to increase the speed of the plots!

Thanks!


texts: Where to download .osm files and how to convert it to .xml?

texts:
I am new to deal with OSM data, and I want to download only .osm files not .pbf. I tried geofabrik.com but it provides .osm.pbf file format

where can I download .osm files?
is there any recommended way to convert .osm files to .xml?


texts: Importing georeferenced 3D objects (Collada DAE) into ArcScene?

texts:
I need to import 3D objects created with Agisoft Photoscan into ArcGis Desktop 10.4.1 (specifically ArcScene). I know how to use the "Import 3D" function to import Collada files as Multipatch files in Arc. There is a problem, however, with the georeferencing of these models as they are in a projected geographic (UTM) coordinate system. This has something to do with the length of the coordinate integers (xxxxxxx.xxxxxx). I know of the "work-around" where one exports two versions of the Collada model from Photoscan - one in local coordinates and one in the geographic coordinate system. You can then import the geographic coordinate model into ArcScene and use the "Replace with Model" function to import the local coordinate version model and fix the bizarre blocky visualization. However, when I do this I notice that there is a slight shift in the location of the models. They are off by a few cm in both x and y, and slightly rotated as well. This means that the models are not correctly positioned in relation to other raster and vector files in the GIS.

I've figured out that I can get around this issue by shortening the length of the UTM coordinates before the decimal when exporting from Photoscan (using the "shift" input). Models with geographic coordinates such as xxx.xxxxxx do not suffer from the same location shift, but this is unfortunate as I have to manually shift every vector and raster in the GIS project to match these shortened coordinates.

Does anyone know of a way to get ArcGIS Desktop to properly display and located full UTM coordinates with imported Collada files?


texts:
I have made a VRT out of a few different single band rasters. I've played around with it in QGIS and it's great and I can do a save as and it'll become a real raster. To do further processing in the pipeline, the files need to be tif files, not a VRT.

How can I batch convert VRTs to real rasters?

Ideally It'd be using command line gdal utilities but I can do it in python too if that's easier.

Or conversely, is there any easier way to just stack rasters (with potentially non-overlapping cell sizes)?


texts: Removing lines from a polygon that has been reprojected to cross the antimeridian

texts:
I have a polygon file in qgis3.0 representing global tectonic plates. I have reprojected the original file in WGS84 (EPSG:4326) to PDC Mercator (EPSG:3832) using only the features that are needed.

Because each feature in the polygon crosses the antimeridian, some (selected) have lines cutting through them:

In the original file the two largest of these features have just one attribute entry, while the other has two but with identical values:

There are probably a number of ways to remove these lines. So far I have tried Simplify and Polygons to Lines to manipulate the nodes but this just causes the outline to become deformed.

Screenshot of invalid geometries:


texts: Does progress bar of ArcPy show when code run from Python window?

texts:
I am learning the usage of SetProgressor. I used the sample code from ArcGIS guide to test, as shown below.

import os
import arcpy

# Allow overwriting of output
arcpy.env.overwriteOutput = True

# Set current workspace
arcpy.env.workspace = r'E:\OilWells Model 20161010\test3\EachRoute'

# Get a list of shapefiles in folder
fcs = arcpy.ListFeatureClasses()

# Find the total count of shapefiles in list
fc_count = len(fcs)

# Set the progressor
arcpy.SetProgressor("step", "Copying shapefiles to geodatabase...",
                    0, fc_count, 1)

# Create a file gdb to contain new feature classes
arcpy.CreateFileGDB_management(arcpy.env.workspace, "fgdb.gdb")

# For each shapefile, copy to a file geodatabase
for shp in fcs:
    # Trim the '.shp' extension
    fc = os.path.splitext(shp)[0]

    # Update the progressor label for current shapefile
    arcpy.SetProgressorLabel("Loading {0}...".format(shp))

    # Copy the data
    arcpy.CopyFeatures_management(shp, os.path.join("fgdb.gdb", fc))

    # Update the progressor position
    arcpy.SetProgressorPosition()

arcpy.ResetProgressor()

When I run the above code in the built-in python shell of ArcGIS, I didn't find the progress bar. Can anyone tell me where is the progress bar showing when run the code?

I am using ArcGIS 10.3.1 Advanced License.


texts: Using proj4js and Leaflet to transform coordinates?

texts:
I'm trying to use proj4js with Leaflet to convert a mouse click event from WGS84 coordinates into British National Grid but when I do this I get NaN for the converted coordinates.

Here's the code:

        function onMapClick(e) {

        var bngprojection = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs";
        var bngcoords = proj4(bngprojection, [e.latlng]);
        popup
            .setLatLng(e.latlng)
            .setContent("You clicked the map at " + e.latlng.toString() + ", or in bng... " + bngcoords.toString())
            .openOn(map);
    }

The BNG projection info I got from: http://spatialreference.org/ref/epsg/27700/proj4/. I can't seem to find much documentation at https://github.com/proj4js/proj4js


texts: Selecting polygons that meet minimum containment of certain raster cell value using ArcGIS Desktop

texts: I have a census block shapefile with a land cover TIFF underneath it. The Raster layer has three different cell values. I am trying to select the census blocks that have a minimum of 50% of cells withing its area with a value of 3. I have access to all extension within ArcGIS Desktop 10.7.


texts: QGIS how to select a set of features in field calculator expression

texts:
I try to connect one feature from a polygon layer to many in a point layer, using geometry generator and expression calculator in the polygon layer.
If I can connect 1 polygon to all the points,
shortest_line(geometry(get_feature('mylayer', 'myattribute', 'myvalue')),
geometry($currentfeature))

I can't connect one polygon to a set of points i.e. I can't select some points with a specific value in a field.
shortest_line(geometry(get_feature('mylayer', 'myattribute', 'myvalue')),
geometry(*only_the_polygons_with_a_specific_attribute_value*))

I tried get_featureand attribute but I failed.
Is there a way to do that ?


texts: Rule-based symbology based on given table

texts:
I would like to ask if it's somehow possible to create script with rule-based symbology in QGIS which would be based on expressions given from separate table (for example CSV, or Excel sheet). I used code given from user eurojam:

def rule_based_style(layer, symbol, label, expression, color):
    root_rule = renderer.rootRule()
    rule = root_rule.children()[0].clone()
    rule.setLabel(label)
    rule.setFilterExpression(expression)
    rule.symbol().setColor(QColor(color))
    root_rule.appendChild(rule)

layer = iface.activeLayer()
symbol = QgsSymbol.defaultSymbol(layer.geometryType())
renderer = QgsRuleBasedRenderer(symbol)
f1 = open("d:/test.csv", "r", encoding='utf-8', errors='ignore')
count = 1
for line in f1.readlines():
   rule_based_style(layer, symbol, str(count), line, "red")
   count = count + 1

layer.setRenderer(renderer)
layer.triggerRepaint()
iface.layerTreeView().refreshLayerSymbology(layer.id())

Problem is that my expression which is for example: "Benzinka" = 'ano' AND "typ_silnic" = 'trunk' AND "Corine" = '121' now gets another quotation marks and the result is: """Benzinka"" = 'ano' AND ""typ_silnic"" = 'trunk' AND ""Corine"" = '121'"


texts: Creating a "VV","VH","VV-VH" composite for each S1 image in loop using GEE

texts:
I would like to create a "VV","VH","VV-VH" composite for each feature in my s1 collection. I am not familiar with loops in GEE.

I filtered the data.

var myCollection = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
    .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
    .filterDate('2015-01-01', '2018-12-31')
    .filterBounds(roi);
print(myCollection);


texts:
I am mapping mountainous areas in QGIS 2.4 and would like to recreate the beautiful blended colour relief found, for example, in the map shown at http://maps.kartoza.com/Boosmansbos/.

While I am able to copy the RGB values of the colours from the key at the side of the map, and map them into my own elevation bands, I am unable to recreate the 'blended' effect between the elevation bands as QGIS produces a 'hard' line between elevation levels.

I'm using the basic mode I suppose which is the menu at the top, Raster - Terrain Analysis - Relief, then mapping the colour bands that way. I thought I might be able to get a 'blend' by overlapping the elevations slightly in each band, but it doesn't seem to let me do that.


texts: grid_canopy failed to create vrt with catalog

texts:
I use this command but it failed to create the vrt file :

chm <- grid_canopy(ctg, 1, pitfree(c(0,2,5,10,15), c(0,1), subcircle = 0.2))

it is possible to add attribute to non create the vrt file to the command grid_canopy ?

the error message :

I know that the question is very short, but i can't send to you the entire las files. This is the error message :

Error in system(cmd, intern = TRUE) :
  'CreateProcess' failed to run 'C:\PROGRA~1\QGIS3~1.10\bin\GDALBU~1.EXE -vrtnodata "-Inf" "C:/Users/po76844/Documents/lidar/la_fresse_la_joux/normalized/grid_metrics.vrt" "C:/Users/po76844/Documents/lidar/la_fresse_la_joux/normalized/14585_Joux-Fresse_14585_Joux-Fresse_Semis_classeL93-IGN69_000001_L93_IGN69_normalized_normalized.tif" "C:/Users/po76844/Documents/lidar/la_fresse_la_joux/normalized/14585_Joux-Fresse_14585_Joux-Fresse_Semis_classeL93-IGN69_000002_L93_IGN69_normalized_normalized.tif"


texts: Can I overlay a population raster with fishnet shapefile?

texts:
I used the dasymetric toolbox to proportion the population of a city based on areas of habitable land. I have the output in a raster file with a corresponding attribute table. My output has a column that has the new population proportioned to each raster zone (there were 42 in total). I have also created a fishnet shapefile.

I'm trying to overlay the raster with the fishnet so I can find out how many people live within each fishnet cell. Though, I am wondering if this is possible in any way.

The software I'm using is ArcGIS 10.4

I am attaching two screenshots. The first is of my raster file with its attribute table and the second is of my shapefile with its attribute table.

My methodology is as follows:
I used the nlcd landcover classification raster and clipped it to my area of interest.I then reclassified the raster where I used 0 for areas I consider uninhabitable by humans, i.e. water, ice, etc. then I used 5, 15, 30 and 65 based on population density. Next, I used the dasymetric toolbox to calculate a new population field (NEW_POP in screenshot 1). I joined the raster file with the dasymetric tool output then I deleted areas of the raster with values = 0. The remaining cells can be seen in Screenshot 1.

After that, I converted my landcover raster to a polygon and did a union of the landcover polygon file and fishnet. Then I deleted fields with values = 0.

The output is seen in Screenshot 2.

Is there a way to find the population (derived from the raster file in the column NEW_POP) within each polygon?

][3]


texts: Is re-projection needed from SRID 4326 (WGS 84) to SRID 4269 (NAD 83)?

texts:
I have US census data against state boundaries stored as  SRID 4269, in an MSSQL DB.

The input data used for calculations against that dataset is stored as SRID 4326.

As far as I can infer from the SpatialReference.org, SRID 4269 is just a subset of SRID 4326 and no coordinate re-projection is necessary.

Do I actually need to re-project from SRID 4326 to SRID 4269, provided that the input data is always within the bounds of the SRID 4269?


texts: Pathfinder Post-processing incorrectly- Data severely shifted

texts: I'm new to Pathfinder (just transferred our license over to my computer). I post-processed my data and exported to ESRI shapefile. My data from the east end of Long Island, New York displays off the coast of Africa. I've tried post processing a few different files, ones that post-processed fine on my coworkers computer, and all have the same issue. I can't seem to find any settings or help documents that address the issue. Any suggestions?


texts: No plugin icon in the toolbox in QGis 2.4

texts:
When i'm loading my plugin into QGis 2.4 the icon doesn't appear in the toolbar. There is only an empty button. I have the icon.png in all the necessary files:

resources.py:

<RCC>
     <qresource prefix="/plugins/colorTransparency" >
        <file>icon.png</file>
     </qresource>
</RCC>

which is compiled with pyrcc4 -o resources.py resources.qrc

metadata.txt:

icon=icon.png

colorTransparency:

 def initGui(self):
    # Create action that will start plugin configuration
    self.action = QAction(
        QIcon(":/plugins/colorTransparency/icon.png"),
        u"Transparency and Random Color", self.iface.mainWindow())
    #connect the action to the run method
    self.action.triggered.connect(self.run)

    # Add toolbar button and menu item
    self.iface.addToolBarIcon(self.action)
    self.iface.addPluginToMenu(u"&TransparencyRandomColor", self.action)

def unload(self):
    # Remove the plugin menu item and icon
    self.iface.removePluginMenu(u"&TransparencyRandomColor", self.action)
    self.iface.removeToolBarIcon(self.action)

The icon.png is in the same directory and has a size of 25px x 25px.

In the plugin installer the icon is displayed.
I also tried to delete it and clean everything but it doesn't work.


texts: What is this Map Algebra expression performing in ArcGIS ModelBuilder?

texts:
I was given a model my someone else, and am trying to decipher exactly what they are doing.

There is a step in the model which has the following expression:

SetNull("%Raster1%"<0,"%Raster 1%")

Is it taking all values of 0 and setting them to null?

What is happening in the second portion of the condition?


texts: I am looking at updating a number of large shapefiles using ogrinfo and UPDATE. I would like to perform a CASE function to update specific fields with specific values. Is it possible to perform CASE functions within ogrinfo UPDATE?


texts: Why do some ArcToolbox tools have red x over them?

texts:
Why is it when you right-click to open some of the tools on the ArcToolbox, such as the TIN to Raster Tool for example, an image of a red x immediately appears over the tool icon?  The tool dialog then never opens.  What am I not understanding?

I'm using ArcGIS 10.2 for Desktop.


texts: Splitting existing edge in postgis topology?

texts:
I try to use topology.ST_NewEdgessplit to  split an existing edge along this edge.
I use the following statement:

  select topology.ST_NewEdgessplit( 'mytopo' , 3 ,(Select ST_Line_Interpolate_Point((select geom from mytopo.edge_data where edge_id = 3), 0.5)));

My expectation is that the interpolated point must be on edge 3 but I get this error message.

  ERROR: SQL/MM Spatial exception - point not on edge
  SQL state: XX000

Why might this go wrong?


texts: Line bearing for wind direction for plume of smoke QGIS

texts:
I'm trying to draw a line representing direction of a plume of smoke for an emergency planning scenario in QGIS. The line needs to be 5km in length and have a bearing of 85 degrees. The closest I have found to achieve this is the CAD like drawing plug-in but this uses angles rather than compass bearings, so 180 degrees would be pointing west rather than south and vice verse.

Is there a tool to draw lines using compass angles and ideally length?

Thanks.


texts: Unable to install QGIS, libnetcdf.so.7 not found?

texts:
I recently updated my system from Ubuntu 14.04 (trusty tahr) to Ubuntu 16.04.4 LTS (Xenial Xerus).

While installing QGIS I get an error message:

/usr/lib/qgis/crssync: error while loading shared libraries: libnetcdf.so.7: cannot open shared object file: No such file or directory

I updated the repositories:

sudo apt-get update

and installed the latest version

sudo apt-get install libnetcdf-dev

While checking the installation of the library I found that libnetcdf.so.11 is installed

Checking out the available repositories for Xenial for ibnetcdf at:
https://packages.ubuntu.com/xenial/libdevel/libnetcdf-dev

I found that the available library is: libnetcdf11 (= 1:4.4.0-2) which corresponds with the libnetcdf.so.11 file I found.

What do I do in order to install libnetcdf.so.7 and fulfil the dependencies for QGIS?


texts:
Is there a way to hide the ArcMap Main Menu either by customizing ArcMap manually or programmatically? I’ve tried to use the customization dialog and disable the Main Menu which didn’t work.

Trying to do it programmatically in the IApplicationStatusEvents_Event.Initialized didn’t work, too.

ESRI.ArcGIS.Framework.ICommandBars commandBars = ArcMap.Application.Document.CommandBars;
if (commandBars != null)
{
    commandBars.HideAllToolbars();
}

Additionally I would like to be able to restore the Main Menu easyly if a support task makes it necessary.


texts:
I am having a really hard time digitizing roads in a map for an assignment in my GIS course. I follow the steps, I open editor, snapping, I open a window to create a new feature. But when I start digitizing, the lines I draw keep disappearing after I double click on them or click finish sketch. I use the auto complete polygon construction tool because I am not making an enclosed polygon I am just making lines along roads.

Does anyone know what issues might causing my lines to disappear and how I can fix them?


texts: Why can't I reproject geopackage layer from WGS84 to ETRS89/UTM (QGIS3.6)

texts:
I accidentally digitized a layer (geopackage) in WGS84 (EPSG 4326). Now I need to reproject it into ETRS89/UTM zone 32N (EPSG 25832). I tried it in different ways which usually work without a problem (saving layer as esri-shapefile or as geopackage layer with the updated projection). But in this case, the resulting layers cannot be displayed (see screenshot with the resulting layer's properties).

What am I missing here?


texts:
We have built an OSM Tile Server following the instructions at

http://switch2osm.org/serving-tiles/building-a-tile-server-from-packages/

The tile server is rendering tiles, but the street names on the tiles are
not readable on high DPI devices, e.g. Nexus 6 at 2560×1440 493dpi.

We would like to modify the tile server so it generates tiles based upon the
client screen's DPI. For example changing the request URL to include a
DPI parameter, e.g.

http://mytileserver.com/osm/client-dpi/12/1050/1522.png

Please let us know how this might be accomplished?


texts: Update raster color ramp with new min/max values using python?

texts:
For a script, I want to update the exsiting min and max of a colour ramp with a new one. In the image you can see what I want to update.

To do so I tried to follow this thread. I tested the following code to update the raster:

min = 1
max = 5

l.renderer().setClassificationMin(min)
l.renderer().setClassificationMax(max)

I get the following error:

AttributeError: 'QgsSingleBandGrayRenderer' object has no attribute 'setClassificationMin'

I can't find any information about this code on the internet. The answer of the thread gives a solution but you HAVE to update the colour as well.

I'm wondering if somebody know how to update the min and max without updating the colour.


texts:
I'm trying to perform an interpolation on SAGA in QGIS but every time I ran the tool (or any other SAGA tool for that matter) I got the following message:

"The following layers were not correctly generated.C:/Users/Ronen/AppData/Local/Temp/processing_6af47249bd814d3e8868defb8ceab637/06a71a2ec1914aa69d2b5a600f18cdef/TABLE.shpC:/Users/Ronen/AppData/Local/Temp/processing_6af47249bd814d3e8868defb8ceab637/3d7ecad0e0fb468a8d6ba4fbdd06b380/SEGMENTS.sdatC:/Users/Ronen/AppData/Local/Temp/processing_6af47249bd814d3e8868defb8ceab637/94135e2aabc549f9b4b2885b81c51cc9/SIMILARITY.sdatYou can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm."

When I check the Log Messages Panel, it displays the following error:

2019-01-03T17:01:35 WARNING warning:C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python\qgis\utils.py:672: DeprecationWarning: the imp module is deprecated in favour of importlib; see the module's documentation for alternative uses mod = _builtin_import(name, globals, locals, fromlist, level)

Does anyone know what that means?


texts: Joining CSV file to Shp file using ogr2ogr and python?

texts:
I am having problems with the code below i.e. parsing the SQL string in python. I keep getting an error. I am joining a CSV file called lookup_scales.csv to a shp file called RC_CATALOGUE_GRANTED. I am running the ogr2ogr command within python and using os.system to call it. It doesn't like the 'lookup_scales.csv' part i.e the apostrophes around it. BTW I was able to run the command successfully within a DOS window so not sure why is fails in the os.system call?

Also the shp file and csv file actually sit in different directories and would like to specify those locations in the ogr2ogr command. At the moment to get the command to work in DOS I have had to put the shp file and csv file in the directory I am running the python code I have.

Code:

os.system('ogr2ogr -sql "select RC_CATALOGUE_GRANTED.*, lookup_scales.* from RC_CATALOGUE_GRANTED left join 'lookup_scales.csv'.lookup_scales on RC_CATALOGUE_GRANTED.SCALE = lookup_scales.Scale" shape_join.shp RC_CATALOGUE_GRANTED.shp')

Above based on the Join CSV file to shapefile thread on here


texts: What types of data can be represented in dot density map?

texts:
I was creating a dot density map for GDP data of African countries. When I submitted the map I got feedback as "this map is incorrect in terms of visualization of data being used"

Each dots means 100 US$. Counting the dots on the map we can calculate the tentative GDP of the country easily. Dots are denser in countries with high GDP where as they are sparse in countries with low GDP. I am wondering how is the map incorrect in term of visualization of data being used?


texts: How to query DEM elevation values in a web application based on cursor/click point location?

texts:
I work for a city government that has LiDAR derived DEM rasters which we would like to make accessible to other departments in the city.  The problem is that many departments don't have anyone with GIS skills and/or licenses.  We are looking to develop a simple, easy to use application that will just return elevation values from the DEM based on cursor click/point location.  We would like to do this without having to serve out the high resolution DEM dataset for speed and ease of use of the application.

We currently have GeoPDFs that can return elevation data as long as users have the TerraGo toolbar installed.  This is expensive and extremely cost-prohibitive so we are looking for a cheap option that can be used by employees and the public alike. We are just looking to query based on XY location of a point which will return the elevation value based off of the DEM.  Has anyone done this or have advice on the best way to proceed?

We do have access to Google Maps Engine (the Elevation API is SRTM so resolution isn't good enough for our purposes), as well as ArcGIS Online.

Also worth mention: We are currently running Oracle Spatial but are transferring to SQL Server 2012 with ESRI's LGIM for our spatial data(geodatabases and feature classes).


texts: How to run multiple tool from different workspace environments

texts:
I've created 1 script to perform 2 specific tasks but it not working correctly due to how the arcpy.env.workspace is set.

The first part of the script deletes and recreates the sub folders in "C:/EIS/SHP Files"

arcpy.env.workspace = "C:/EIS/SHP Files"

The second part of the script looks inside a particular sub folder in "C:/EIS/SHP Files/Master folder" and copies Pipe data back into the newly created folders

arcpy.env.workspace = "C:/EIS/SHP Files/Master Folder"

Problem: I can't reference 2 environments at once in a script. The script will give out an error message since it can't look in 2 paths at once. So how can we do this?

import sys
import platform
import imp
print "Importing arcpy... this may take a moment\n"
import arcpy
from arcpy import env

# Set workspace
arcpy.env.workspace = "C:/EIS/SHP Files"

# Creates the county folders
out_folder_path = "C:/EIS/SHP Files/"

# Set local variables
# List of folder values in load order
folderList = ["1","2","3","4","5"]

# Store all the folder names in a dictionary with the folder code as the
# "key" and the folder description as the "value" (folderDict[code])
folderDict ={"1":"OH_ADAMS","2":"OH_ALLEN","3":"OH_ASHLAND","4":"OH_ASHTABULA","5":"OH_ATHENS"}

# Wayne Pipe
in_data_WayneMain =  "WayneMain.shp"
out_Geauga_WayneMain = "C:/EIS/SHP Files/OH_GEAUGA/WayneMain.shp"

# The Code below will utilize the [arcpy.CreateFolder_management] tool in
# order to rebuild the directory for C:\EIS\SHP Files\.. County Folders

# Process: Add valid material types to the domain
# use a for loop to cycle through all the domain codes in the dictionary
for code in folderList:
    if arcpy.Exists(out_folder_path):
        arcpy.Delete_management(folderDict[code])
        arcpy.CreateFolder_management(out_folder_path, folderDict[code])

# The Code below will utilize the [arcpy.Copy_management] tool in order to
# copy files from the C:\EIS\SHP Files\Master Folder to C:\EIS\SHP Files..
# County Folders

# Wayne Pipe (Copied) To => County Folders
arcpy.Copy_management(in_data_WayneMain, out_Geauga_WayneMain)


texts: What is the best raster format for large NoData areas?

texts:
Is there a good raster format for dealing with a large number of NoData cells?

I do a lot of flood modeling and mapping, and recently I've been using rasters to show predicted flood depth and/or elevation for large areas, often for an entire county.

I've just been using the default ESRI binary grid format, but I'm wondering if there might be something more appropriate. These rasters use a huge amount of space (especially at cell sizes of 5-10 ft), and the vast majority of the cells are NoData, so they're just going to waste.


texts: How can I change the Datastore of a Published Layer in Geoserver?

texts:
I have Published a Layer on Geoserver using Shapefile as the datastore. This data has now been stored in a PostGIS database as a Spatial Table. This PostGIS database has been added a new DataSource.

Is it possible to just change the DataStore of a published layer, so that it will now point to this new DataSource? This will avoid having to set all the setting once again.

How can I do this?


texts: Dissolving points based on easting and northing with statistic type first as aggregate function using PostGIS?

texts:
I have two points features having point_ID, Easting and Northing fields and want to dissolve based on Easting and Northing with statistic type aggregate as "first" for point_ID, as you do in ArcMap.

How can I achieve this using PostGIS?


texts: FME transform query result to JSON and then XML

texts:
I execute query on my postgresdb to get result as a json: description is here http://rextester.com/ZULI42144

It looks nice but now I would like to convert it to xml. I decided to use fme to query postgres just like in example and now I have noe idea how to convert the result to xml and write it to file.

Also talking about json which I also need (I need result in JSON and XML format) what I achieved is that I simply executed postgres query end then wrote results using simple text writer so I suppose that fme doesnt really know that the query result is a JSON :)


texts: Google Earth Engine: trying to find areas of multiple water patches in an image

texts:
I am trying to measure the individual areas of many patches of water in a single image, rather than their total area. The connectedPixelCount function won't work for me because of its upper limit of 1024 pixels for the regions. However, using the reduceToVectors function seems to works to identify contiguous regions of water.
The issue is that I am unable to calculate the areas of the resulting polygons.

This part of the code appears to work OK, setting the basin boundary,
getting the image to be used, finding the water regions and plotting them.

// get basin boundary
var fc =
ee.FeatureCollection('ft:v:1L6-BUK8GNDUU-335Ju8D9Gf9mXiOuPKigGSFcf0w');
var startDate = ee.Date.fromYMD(2015, 5, 1);
var endDate = ee.Date.fromYMD(2015, 5, 30);

// set region to be viewed
Map.setCenter(-106.086940, 52.207365, 13);   // St Denis

// get imagecollection, filter by start and end date
var monthly=ee.ImageCollection('JRC/GSW1_0/MonthlyHistory')
  .filterBounds(fc)
  .filterDate(startDate, endDate);

// plot basin boundary
Map.addLayer({
   eeObject: fc,
   name: 'St Denis Basin',
});

// get single image, clip to basin boundary
var current = ee.Image(monthly.first())
 .clipToCollection(fc);

// get water, defined as value 2
var water = current.eq(2);
water = water.updateMask(water.neq(0));

// get water regions as vectors
var vectors = water.addBands(current).reduceToVectors({
  geometry: fc,
  scale: 20,
  crs: current.projection(),
  labelProperty: 'water',
  geometryType: 'polygon',
  eightConnected: true,
  reducer: ee.Reducer.mean(),
});

// plot water vectors
var empty = ee.Image().byte();
var fills = empty.paint({
   featureCollection: vectors,
});
Map.addLayer(fills, {palette: '0000FF'}, 'colored fills');

This section does NOT work:

var areas = function(feature) {
// Compute area from the geometry.
  var area = feature.geometry().area();
  return feature.set('area', area);
};

// Map the difference function over the collection.
var allAreas = ee.Number(vectors.map(areas))

// Print the result.
print(allAreas);

When I run the code, I get the error:

Number (Error)
Error in map(ID=+411164+143465):
Geometry.area: Unable to perform this geometry operation. Please specify a non-zero error margin.

The ID refers to the first polygon

: Feature +411164+143465 (Polygon, 2 properties)
type: Feature
id: +411164+143465
geometry: Polygon, 9 vertices
type: Polygon
coordinates: List (1 element)
0: List (9 elements)
  0: [-106.13101008465881,52.22454718910092]
  1: [-106.12885412797692,52.22454718910092]
  2: [-106.12885412797692,52.22490651521456]
  3: [-106.1293931171474,52.22490651521456]
  4: [-106.1293931171474,52.22508617827138]
  5: [-106.13029143243152,52.22508617827138]
  6: [-106.13029143243152,52.22490651521456]
  7: [-106.13101008465881,52.22490651521456]
  8: [-106.13101008465881,52.22454718910092]
  geodesic: false
  properties: Object (2 properties)
    mean: 2
    water: 1


texts: Displaying specific version of geodatabase on AxImapControl using ArcObjects in C#?

texts:
I am working on a sample ArcGIS Engine project which tries to map two versions difference on an AxImapControl

i created two mapcontrol which first one loads the MXD file and for  second one i make a deep cloning copy of map object using IObjectCopy interface and then change the version of copied object and the load the map in mapcontrol.

but when i load map to second mapcontrol, two mapcontrols are showing the target version which i want to display on second map.

deep cloning in second map's parent form constructor

   IObjectCopy esriObjectCopy = new ObjectCopyClass();

        System.Object copyTarget = axMapControl1.ActiveView.FocusMap;

        esriObjectCopy.Overwrite(pMap, ref copyTarget);

change version after deep cloning

 ESRI.ArcGIS.Carto.IBasicMap basicMap = axMapControl1.ActiveView.FocusMap as ESRI.ArcGIS.Carto.IBasicMap;
        IChangeDatabaseVersion oIChangeDatabaseVersion = new ChangeDatabaseVersionClass();
        ISet set = oIChangeDatabaseVersion.Execute(_ParentVersion, _ChildVersion, basicMap);


texts: Assigning new QgsRasterBlock object to QgsRasterDataProvider object?

texts:
Next code, read and print original and modified raster pixels values by using a QgsRasterBlock object ('block' method of QgsRasterDataProvider). The objective is to obtain a new raster with the new block values. However, these new values are not overwritten in the provider; as it can be observed at next image (with the help of Value Tool plugin).

Complete code is:

layer = iface.activeLayer()

provider = layer.dataProvider()

extent = provider.extent()

rows = layer.height()
cols = layer.width()

xsize = layer.rasterUnitsPerPixelX()
ysize = layer.rasterUnitsPerPixelY()

print rows, cols

block = provider.block(1, extent, cols, rows)

print "original block"

for i in range(rows):
    for j in range(cols):
        print block.value(i,j)

for i in range(rows):
    for j in range(cols):
        if block.value(i,j) == 4:
            block.setValue(i, j, 25)

print "modified block"

for i in range(rows):
    for j in range(cols):
        print block.value(i,j)

pipe = QgsRasterPipe()

pipe.set(provider.clone())

rasterWriter = QgsRasterFileWriter("/home/zeito/pyqgis_data/test_write.tif")

xSize = provider.xSize()
ySize = provider.ySize()

CRS = layer.crs()

error = rasterWriter.writeRaster(pipe, xSize, ySize, provider.extent(), CRS)

if error == QgsRasterFileWriter.NoError:
    print "Raster was saved successfully!"

else:
    print "Raster was not saved!"

After running the code at the Python Console of QGIS, it can be corroborated that block object was modified as expected.

3 3
original block
10.0
2.0
2.0
10.0
3.0
2.0
8.0
4.0
4.0
modified block
10.0
2.0
2.0
10.0
3.0
2.0
8.0
25.0
25.0
Raster was saved successfully!

When I try to create a new instance of QgsRasterDataProvider, I got that "cannot be instantiated".

TypeError: qgis._core.QgsRasterDataProvider represents a C++ abstract class
and cannot be instantiated

How can I overwrite correctly the new block values in the original (or a new) provider?


texts: Getting list of QgsPoint from polygon Layer in PyQGIS

texts:
How do I manage the points in a layer?

from qgis.core import *

def run_script(iface):
    layer = QgsVectorLayer(layerType.Polygon+'?'+crs.psg_4326,
                       'la-florida',
                       vectorLayerType.memory)

In this moment I need to get all vertex from the polygon which represent the region bound, but I have not find the documentation about how to get it.


texts: Why are the 500 meters not 500 meters? And how do I add 500 meters in x/y direction to a lat/lon point?

texts:
I wanted to create a square cell on a map with size 500m times 500m. Therefore, I wrote the following code:
from pyproj import Transformer

bl_lat, bl_lon = 40.64423086189233, -74.0307935018417
trans_coord_to_meter = Transformer.from_crs('epsg:4326', 'epsg:3857')
bl_t = trans_coord_to_meter.transform(bl_lat, bl_lon)

print(trans_coord_to_meter.transform(bl_t[0], bl_t[1], direction='INVERSE'))  # blue
# 40.64423086189233, -74.0307935018417 (lat, lon)
print(trans_coord_to_meter.transform(bl_t[0]+500, bl_t[1], direction='INVERSE'))  # red, + x_axis
# 40.64423086189233, -74.02630192542111
print(trans_coord_to_meter.transform(bl_t[0], bl_t[1]-500, direction='INVERSE'))  # green, - y_axis
# 40.640822707278815, -74.0307935018417
print(trans_coord_to_meter.transform(bl_t[0]+500, bl_t[1]-500, direction='INVERSE'))  # yellow, + x and y axis

I thought that this adds 500m in both dimensions (maybe I misunderstood the function). But, when I check the result, I observed the following:

So, how do I really get points an actual distance of 500m in between? (for x and y)
I want to keep using Python, but I'm flexible with the library.


texts: How to apply a radial fill to a voronoi polygon centred on its original home point in QGIS?

texts:
I'm working with a point layer (airports from OpenFlights), and I want to create a map of distance to the nearest airport.

I thought I might be able to do this by using radial fill on the voronoi polygons.

I created the voronoi polygons. These included the original fields from the airports layer, so I have (lat, lon) of the airport for each polygon.

When it comes to styling, I tried to create a radial fill centred on the airport. It's very easy to center on the centroid, but the centroid is nowhere near the airport in most cases.

I set up a radial fill with reference points (0.5,0.5) and (1.0,1.0), mode: reflect.

You can see here the result.. how can I centre it so the brightest point in the radial fill for each polygon aligns with its airport?

There doesn't seem to be an option to set a centre point coordinate. The data-driven option to the right of the centroid isn't for setting up a coordinate, it's a binary field for the centroid setting.

I could try computing an offset (between the centroid and the airport) but there aren't any data-driven options for the offset fields, so they would need to be set manually.

Shapeburst fill isn't suitable for this (although it does give some interesting results!)

I've also looked at raster and svg fill, but again these appear to offer no way to set the centre coordinate (I suspect they align using the bounding box)

Is there some way to do this with QGIS? Am I overlooking something?
Also, is there way to prevent QGIS from scaling the gradients to fit the bounding box? I've tried switching from 'object' to 'viewport' but the polygons all vanish.

Dmitry helped fix the first question - aligning up the gradients to the airport locations

I added 6 new real fields to my polygons with the field calculator...

xmin -> x_min(bounds($geometry))
ymin -> y_min(bounds($geometry))
xmax -> x_max(bounds($geometry))
ymax -> y_max(bounds($geometry))
xoff -> (airportx - xmin)/(xmax-xmin)
yoff -> 1 - ((airporty - ymin)/(ymax-ymin))

I then used the values of xoff and yoff in the data-driven x1 and y1 fields. The bounding box appears to be in the bottom-left, not the top-left, hence the "1-" in the last function.

Results now look a lot better...

This still leaves me with the problem of the best option for x2 and y2; for now, I've set the second control point to (1,1). Need to have a think about that.

Given that this is pegged to the bounding box, I suspect the answer to the second question (can I avoid the gradient being scaled) is 'no' :)

So I'll probably do this as a python gdal raster script.. also has the advantage that I can use PostGIS to calculate spheroid distances.


texts: Can I change the default fill style of layers?

texts:
For every layer I add, I have to change the fill style.

Is it possible to define the default fill style on the project level?


texts: Path Distance calculation limited to downhill only

texts:
I am trying to calculate a path between two points with a relative elevation difference of about 100m and 3km separation. The points represent the start and end of a gravity fed waste water pipe. Essentially the path should follow the contours and have a consistent gradual fall.

I have run the Esri Path Distance tool using Slope and DEM as cost and input surface (respectively) with the DEM as the vertical input as well using default values. The resultant path minimises slope but I also want it to only travel downhill.

Can someone explain how to control the vertical factor parameters to only calculate a downhill path. The Esri help has not really described this.


texts: Is it possible to have self intersecting polygons in shapefile?

texts:
GIS vector data for example shapefile contains polygonal data.

Lets assume, we have polygon A which has self-intersection which means that the edges in this polygon intersect with each other.

Lets say if we want to intersect A with polygon B, then the steps of computing the overlaid polygon differs whether A is a polygon with self-intersections or not.

Do the polygons in one map layer have self-intersections or not?

I searched for this question on internet and I found couple of instances where people complain about error caused while importing data due to the presence of self-intersecting polygons.

Also I got the impression that self intersecting polygons are not possible in certain file formats.


texts: Python scripting in arcmap - is there an easy way to re-compile py-file

texts:
in Arcmap 10.2.2 I am editing downloaded toolboxes (not python-toolboxes) by editing the python script in a texteditor. Whenever I restart Arcmap, the Compiled Python File gets updated and I can use the changes I made in my tool in Arcmap.

However, the Compiled Python File does not update without restarting Arcmap (which is time-consuming), why I can not make use of the changes I made without restarting.

Is there an easy way to update / re-compile the  Compiled Python File without restarting Arcmap?


texts:
Is there a way to create a centreline in QGIS from two other lines (such as a river bank) which won't always be equally apart?
I am aware that in ArcGIS Desktop there is Collapse Dual Line to Centerline which sounds like it is what is needed, but I have not found an equivalent in QGIS.


texts: Splitting a map into small grids based on coordinates

texts:
I have an API that requires coordinates as rectangular search box. (x1,y1) is the upper left corner, and (x2,y2) is the lower right corner of the search box.

I am looking into a large area therefore have to split my area into smaller zones, namely grids. I wonder what is the best way to split my entire map into say, 1000 number of smaller grids and feed each smaller grid corresponding coordinates to the API. In other words, get, $x_1_n$, $y_1_n$ and $x_2_n$ and $y_2_n$ where n is the grid id?

I am new in the field, but understand these calculations are on the basis of a spherical Earth. Can you provide any hint or even keywords that will get me into a right path and eventually write my script?


texts:
The Encoding and the Schema tag for ComplexData in a DescribeProcess-Response document is an optional parameter (defined in the WPS Specification).
In GeoServer, the Response document doesn't contain these parameters.
Also, for a literal input parameter, the parameter "DataType" isn't shown.

What do I have to change in GeoServer, that these parameters are shown in the response document? (Respectively, how can I edit the parameters?)

Thanks for helping.


texts: XML Document points to wrong location after using arcpy.BatchProject()

texts:
I'm trying to reproject a number of shapefiles representing contours from NAD_83_Albers to NAD_83_UTM_Zone_9N using arcpy.BatchProject(), but I have to use python/regex to preserve a particular file architecture.

The script works for all practical purposes, however one thing concerns me. When the script finishes, the XML file for the new shapefile points to the location of the original file as the source and lists the old projections instead of the new ones.  The script appears to have simply copied the contents of the old XML file into the new XML instead of showing the proper values.  When I go into the properties of the shapefile in Arc it says the projection is UTM rather than Albers as it should be, but the XML file does not line up with this.

Any ideas?

Here's my code.

import arcpy, os, re
rcont = re.compile(".*\contours")
inputDir = arcpy.GetParameterAsText(0)
outputDir = arcpy.GetParameterAsText(1)
for path, dirs, files in os.walk(inputDir):
  holder = path
  if rcont.match(holder):
    splitList = re.split("[\\\]", holder)
    cell = splitList[4]
    trimName = splitList[5]
    featureType = splitList[6]
    newDir = outputDir + "\\" + cell
    newDir2 = newDir + "\\" + trimName
    newDir3 = newDir2 + "\\" + featureType
    if not os.path.exists(newDir):
        os.mkdir(newDir)
    if not os.path.exists(newDir2):
        os.mkdir(newDir2)
    if not os.path.exists(newDir3):
        os.mkdir(newDir3)
    arcpy.env.workspace = holder
    fileList = []
    for files in arcpy.ListFeatureClasses():
        fileList.append(files)

    res = arcpy.BatchProject_management(fileList, newDir3, sr)
    if res.maxSeverity == 0:
        arcpy.AddMessage("projection of all datasets successful")
    else:
        arcpy.AddMessage("failed to project one or more datasets")


texts: How to fix strange rendering behaviour of OSM Data in Mapnik?

texts:
My current setup so far is a running Apache webserver with tilecache.cgi script and a mapnik rendering engine connecting to my postgis database containing my OSM data.

The problem is very easy to describe. All roadnames, roadnumbers, descriptions are to close to each other. I set up tilecache with MetaTile=true, TileSize=8x8 and it doesn't change. What am I doing wrong?

Thanks!


texts: Use attribute of all filtered features of one layer as label in the Print composer

texts:
Is it possible to filter some features in one layer and use those filtered features as labels in the Print composer of another layer?
I am trying to make an atlas. I have two layers, one which has the province data, and another which has the city data. Both these layers have 9-digit codes; the city data has something like 140106000, while the province has something like 140100000. I would like to filter the features in the city layer to match the code as 1401. All features in the city layer with codes that match 1401 I would like to display their "name" attributes in the print composer (so that all filtered city names appear in the print composer, like a list).
I do not know of any way to do this.


texts: Can ArcGIS web appbuilder work offline without Portal

texts: I have used web app builder for ArcGIS Online. Is there a way for it to work offline without Portal? Consuming ArcGIS Server REST API


texts:
I have a PostgreSQL/PostGIS database located on a server.
I have added shapefiles (with the "PostGIS shapefile import/export manager") and raster (with "raster2pgsql"). It works well.
However, when I try to save layer styles in this database (through QGIS) it works well for shapefiles but for raster it is only possible to save the style on local machine in a *.qml or *.sld file (and therefore not in PostgreSQL as it can be done for shapefiles and stored in layer_styles).
Is there a way to save style for rasters in PostgreSQL/PostGIS ?
(I am a rookie with PostgreSQL/PostGIS).


texts: How to obtain elevation info from VRT file?

texts:
I have a VRT mosaic file I made from thousands of ESRI GRID files.

It looks something like this.

<VRTDataset rasterXSize="111831" rasterYSize="94785">
  <SRS dataAxisToSRSAxisMapping="1,2">PROJCS["unnamed",GEOGCS["Unknown datum based upon the GRS 1980 ellipsoid",DATUM["Not_specified_based_on_GRS_1980_ellipsoid",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6019"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4019"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",24.3333333333333],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.999941176470588],PARAMETER["false_easting",656166.666666667],PARAMETER["false_northing",0],UNIT["Foot_US",0.304800609601219],AXIS["Easting",EAST],AXIS["Northing",NORTH]]</SRS>
  <GeoTransform>  6.9303750000000000e+05,  2.5000000000000000e+00,  0.0000000000000000e+00,  9.5992000000000000e+05,  0.0000000000000000e+00, -2.5000000000000000e+00</GeoTransform>
  <VRTRasterBand dataType="Float32" band="1">
    <NoDataValue>-3.4028234663852886e+38</NoDataValue>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">be_rasters/692500_725000</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="785" RasterYSize="1000" DataType="Float32" BlockXSize="256" BlockYSize="16" />
      <SrcRect xOff="0" yOff="0" xSize="785" ySize="1000" />
      <DstRect xOff="0" yOff="92968" xSize="785" ySize="1000" />
      <NODATA>-3.4028234663852886e+38</NODATA>
    </ComplexSource>
  </VRTRasterBand>
</VRTDataset>

Of course, my real file has thousands of ComplexSource elements. When opened in QGIS, it shows the elevations of every property in a U.S. county.

Here's how part of it looks in QGIS.

How do I obtain the elevation levels of each property, in feet or meters, from the VRT mosaic?


texts:
I'm re-opening this question (Open JPEG2000 (Sentinel 2) in R). I'm working on a Mac OS X Yosemite 10.10.5.

I would like to be able to read jp2 format images with R without translating it to .tif format or anything of that kind.

I followed what have been explained in this question, but I'm still not able to open jp2 images.

library(rgdal)
Le chargement a nécessité le package : sp
rgdal: version: 1.2-7, (SVN revision 660)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
Linking to sp version: 1.2-4

Then :

s04 <- readGDAL("/Users/martinmonziols/Documents/ENSAE/Stage TheGreenData/L2A_T30UYA_20170509T105621_B04_10m.jp2")
Error in .local(.Object, ...) :

So, I checked the drivers

library(gdalUtils)
gdalDrivers()

It seems that I don't have any driver for jp2 files. But what is strange to me is that when running these lines in Terminal

MacBook-Air-de-Martin:~ martinmonziols$ gdalinfo --version
GDAL 2.1.3, released 2017/20/01

And

gdalinfo --formats

I have in particular this line in the table of all drivers:

JPEG2000 -raster,vector- (rwv): JPEG-2000 part 1 (ISO/IEC 15444-1), based on Jasper library

which means that I have somewhere a driver for jp2 format. But R is not speaking to this one since the versions differ and this driver does not appear in the list I have in R.

What am I missing here ? how can I force R to use the right gdal version instead of the older one ?


texts: very basic question on extracting data from .tif raster layer in R- projection not in lat/long

texts:
I have data on N deposition in the US, which can be found in a .tif file from a zipped directory, which can be accessed here: nadp.sws.uiuc.edu/maplib/grids/2006/TotalN_dep_2006.zip.

I loaded the data into R using the raster package:

require(raster)
r <- raster('/path/to/dep_totalN_2006.tif')

I then have some coordinates for which I would like to extract values from the raster. Here are 6 sites for testing, in dput format:

test.dat <- structure(list(latitude = c(46.414597, 46.137664, 42.258794,
44.287538, 46.567187, 46.205438), longitude = c(-86.030373, -85.990492,
-85.847991, -85.806588, -87.954285, -87.481934)), .Names = c("latitude",
"longitude"), class = "data.frame", row.names = c(NA, 6L))
> test.dat
  latitude longitude
1 46.41460 -86.03037
2 46.13766 -85.99049
3 42.25879 -85.84799
4 44.28754 -85.80659
5 46.56719 -87.95428
6 46.20544 -87.48193

Generally, extracting data is pretty straightforward from this point on. I would grab data from the layer for these point using something like this:

points         <- cbind(test.dat$longitude,test.dat$latitude)
test.dat$out   <- extract(r, points)

This doesn't work however. It does not error, it just produces a vector of NA values. I can plot the raster layer, which gives a clue to whats going wrong:

plot(r)

Clearly the x/y coordinates of the raster layer are not in latitude/longitude, and I expect this is why my extract command fails to extract any data from the layer for my list of points. Apologies if this is very basic and has been answered elsewhere, but my searching hasn't yielded a simple answer that came across clearly. I'm experienced with R, but very light on my use of spatial data within R. How can I accurately project the raster file, r such that my call to extract generates the appropriate values for my sites within test.dat based on their latitude/longitude?


texts: How to view error feedback in qgis when testing connection to postgres database?

texts:
Using qgis 1.8 on Ubuntu 10.04 with Postgres 8.4.13

I am creating a new (and my first) Postgresql connection in qgis.

I can connect to the user and the db using the pgsql client.

jlam@lusk:~/programmation/Quantum-GIS/build-master/output/bin$ psql -U qgis -W qgis -p 5432
Mot de passe pour l'utilisateur qgis :
psql (8.4.13)
Saisissez « help » pour l'aide.

qgis=> \q
jlam@lusk:~/programmation/Quantum-GIS/build-master/output/bin$

But when testing the same parameters within the PostGIS connection creation dialog box, I get connection failure:

"La connexion échouée - Réglez la connexion et réessayez"

("Connection failed - Fix the connection and try again")

My pg_hba.conf file includes the line:

# "local" is for Unix domain socket connections only
local   all         qgis                          md5

Is there a way I can get qgis to send me back the error code and error string it gets from postgres, if its connecting at all, or to increase debug output on the console?


texts: How to determine relative position of a new waypoint in an array of waypoints?

texts:
I hope this is not a ridiculously stupid question but I can't seem to come up with an immediate algorithm for this:

I have a list of (lat,lng), which represent a list of waypoints on a map. At any point the user can click on the map and create a new waypoint corresponding to that new (lat,lng) coordinate. The waypoints correspond to a flight path to follow, and I potentially have thousands of waypoints already defined. The new waypoint needs to be inserted in between the 2 closest waypoints already defined in the list of waypoints. This needs to be done as efficiently as possible.

One option would be to calculate the great circle distance between the new waypoint and the all the previously defined waypoints, and then add the new waypoint in between the 2 closest.

I was wondering if there is a better algorithm for this? Any pointers to efficient/fast algorithms would be very helpful


texts: QGIS export to shapefile inflates real field size and value

texts:
At 3.8.3, when I Export > Save As a Geopackage point layer to a shapefile, fields with type = double are output with 15(!!!) significant digits.  Here are the before and after field widths:

Input:
Type:    double
Length:   10
Precision:   2

Output:
Type:    double
Length:   23
Precision:  15

Furthermore, the output values are modified by a slight amount.  Here are some examples:

Input    Output
150.62   150.620000000000005
150.51   150.509999999999991
221.86   221.860000000000014

How do I maintain the source length and precision?

EDIT:

Based on the comment kindly provided by user30184, the QGIS shapefile output that my OP described appears - unfortunately - to be the default action.

Now I'm curious about the passage that user30184 quoted, particularly the phrase "...Real (floating point) fields without an explicit width..."  This implies that my fields did not have explicit widths.  My new question is, how do I make my real field widths explicit, either by modifying the source layer, or in the GDAL output process, or some other method?