texts: Setting selection border color in red and save the map of the resulting selection

texts:
I'm attempting to create a plugin in QGIS 3 in order to zoom on a selected feature and get the map picture exported in .jpg format.

I've been able to perform the zoom and selection of my selection, and now I want to add some color border to the selected feature using QgsHightlight :

def highlight_features():
    for i in couche_parcelles.selectedFeatures():
        h = QgsHighlight(canvas,i.geometry(),couche_parcelles)
        h.setColor(QColor(255,0,0))
        h.setWidth(10)

#================================Map Creation==================================
#project = QgsProject.instance()

mapSettings = canvas.mapSettings()
mapSettings.setOutputSize(QSize(1600,1100))
mapSettings.setOutputDpi(200)

render = QgsMapRendererSequentialJob(mapSettings)
render.start()
render.waitForFinished()

dir = r"C:\Users\[...]\Desktop"

img = render.renderedImage()
img.save(dir + "/" + "selection.jpg",'jpg')

It does make a pretty highlight border on QGIS software like this :

But then, once I get the map rendered, here's what I get generated by the script :

My question is, how can I get the highlight extracted with my selection on the map ? Because I'm not sure it has something to do with my canvas and I don't know if it's possible to embed my QgsHighlight in my canvas either

I've already seen this answer : Setting selection color transparent and border color red using PyQGIS?
Which helped me to get the highlight on my selection borders but I can't figure out how to integrate it in my map rendering


texts: Running ArcGIS field calculation using Python with algebraic expression using other field values as variables?

texts:
I'm trying to run a field calculation using multiple fields in a table as well as independent numeric values with python code in the field calculator window.  I've included the first 7 lines of the code below.  I'm not getting any syntax errors, but I'm also not getting any results.  The other fields used as variables include: Appraised_Value, Ag_value, and TaxRate_Municipal  Not sure how to fix it.

Pre-Logic Script Code:

def Reclass( Exemp_Code_Int ) :
 if (Exemp_Code_Int == '0') :
  return (Appraised_Value / 100) * TaxRate_Municipal
 elif (Exemp_Code_Int == '2') :
  return (Ag_Value / 100) * TaxRate_Municipal
 elif (Exemp_Code_Int == '4') :
  return 0

Codeblock looks like this:

 Reclass( !Exemp_Code_Int! )

Screenshot too:


texts: Spatialite return XY coordinate in Double Precision

texts:
I'm trying to return the first vertex of a polyline in QspatiaLite. The spatial function x() is used, but it returns an integer coordinate value. How can I return the coordinate value as double precision? This is my query:

SELECT id, AsText(StartPoint(Geometry)), X(PointN(Geometry, 1))
FROM shp_lane

This is the result:


texts:
When I try to open my .qgs file I get the following error:

It says that the project can't be opened, and that something is wrong when trying to read the project file, this is my translation.

Can I do anything to "fix" the file so I can open it?

No, I have no backup, just a .qgs~, does that help?
I use QGIS 2.12.1 and Windows 10, and the file I'm trying to open is also made in 2.12.1.


texts: Expected result for ZonalStatisticsAsTable MEAN Statistic Type?

texts:
I'm getting different results when using the ZonalStatisticsAsTable tool in 10.4.1 when compared to earlier versions (10.2 and 10.3).

ArcGIS 10.2:

ArcGIS 10.3:

But in ArcGIS 10.4.1, the same data produces this table as a result:

For each test, I used the same statistic type and checked the option to ignore NODATA in the calculation.

Has the earlier versions been incorrect, and 10.4.1 is correct? Or is 10.4.1 incorrect?

I have a set of reproducible test data here:
https://bitbucket.org/CumminsJP/zonalstatisticsastable-issue-104/src


texts: How to change geotiff file using the header?

texts:
I have some geotiff file and I need to change the position: 500 meters to east and 500 m to south. Resolution is 12,5 meters.

I've tried using GeoTIFF Tools GUI, and I've created .gtf file, but I don't know what I should write in.

Also I've tried with GeoTiffExaminer, but after I wrote in X vector X and Y vector Y values nothing changes.


texts: Automate end-snap of lines mid-way between nodes?

texts:
I have a polyline feature class of field-collected (Trimble GeoXT with ArcPad10) stream segments.  For whatever reason (still troubleshooting why), my employee returned from the field and the majority of the endpoints were not snapped (despite occurring within distances of the snapping distance I'd specified on the XT).

To automate generation of a continuous series of stream segments, I already ran the "Snap" tool, but it only moved one endpoint and left the other one fixed.  In this scenario, I cannot assume that one end-point is more accurate than the other, so I would like to average them.

Is there a command (or maybe an environment setting within "Snap") I'm missing that would move endpoints for adjacent lines and snap them half-way along a straightline distance in between?

Alternatively, a command that generates a point midway between ends that I could subsequently snap lines to would probably work as well.

I have an ArcInfo license v10 SP4.

Thanks,
Will


texts: Importing CSV numbers formatted as string into ArcMap?

texts:
I'm importing a .csv that contains data that I am joining to a shapefile. The value that I am joining is shown in the first column of the image below. When I import the data into ArcMap the data in the second column is converted to type 'string' instead of 'integer'.

How do I import the data so that numeric values are of type integer 'instead' of 'string'?


texts: How to get started on web GIS development with Java?

texts:
I am very new for GIS, I know java web development.

I want to develop a GIS web application in Struts, can you suggest me how to start?

My requirement are as -

Any API/Tool should be open source.
It should work on HTTPS.

Any tutorial?


texts: Openlayers using Google Maps as base layer: possible to zoom out farther than max extent?

texts:
See attached picture. This is how it is when the page loads. I'm currently using

map.zoomToMaxExtent();

Everything is working pretty well but I'd like to zoom out more/see the whole world using OSM/Google Map Layers: is this even possible?

Note that the zoom bar on the left is set to the lowest bar when the page loads (if you can't see the map well, just right click the image and paste into a new browser tab).

This post was somewhat helpful but even when I do MIN_ZOOM_LEVEL that doesn't do much.

Also note I am not setting any bounds, extents, or resolutions. Toying with those did not seem to do the trick either.

Any input is appreciated.

Edit: here's a codeshare if you want to take a look: http://codeshare.io/GLC4i

UPDATE: http://i.imgur.com/mSpTchY.png This is not perfected yet (and it makes me think something may be wrong with my OL code itself) but currently I am using this

var extent = new OpenLayers.Bounds(-1.32,81.71,-1.18,-61.80).transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:3857"));

I tried

var extent = new OpenLayers.Bounds(-20037508.34279,-10018754.17139,20037508.34279,10018754.17139);

(both with and without transform) but this did not do anything.

UPDATE 2: http://i.imgur.com/XzuGxrn.png

wrapDateLine: false works fine for now: can zoom out farther, pan around, etc etc.

Now I would like to cut off everything to the left of the date line on the left side of the map (the 'first' portion), and cut off everything to the right of the date line on the right side of the map (the '3rd' portion) WHILE still keeping the map full screen. I've tried bounds, resolutions, but no go. Is this even feasible?

Something like this but it should take up the full screen: possible to do? See below:

http://i.imgur.com/feBEnx3.png


texts: generate a polygon that only keep its outside boundary

texts:
I have 25 different polygons as show in fig. below. These small  polygons can form a bigger polygon without overlaps.

My goal is, according the data of 25 smaller polygons, to generate the bigger polygon that only keep its outside boundary.

Now I have selected the 25 polygons and saveas another file, but the innner lines are also included.

I have tried the layer-property-symbol level method, but there is no such button in my .shp file. I have no idea why.


texts: Joining single column in ArcGIS Desktop with Spatial Join

texts:
I'm trying to automate a rather annoying workflow in ArcGIS 10.4 with a Python script. My goal is to extract data from a polygon shapefile into another polygon shapefile which is arranged as a raster (see below). For that I use the tool intersect. Since I´m just interested in the green area in every square I than calculate the geometry of the intersection polygon. Finally I want the area column of the intersection polygon to be joined to the the "grey" polygon file and for that I would use spatial join. All together this works so far but I want just the area column to be joined and spatial join is joining all columns. I can manually delete these but since I'm having multiple files I don't want to do this.

With the help of Google I found fieldmapping but this is not really working and in my opinion its too complicated for what I want.

Is there a simpler approach maybe more Python-related?

It would be great if I could just join only the needed column in the first place.

Has anyone an idea on how I can solve my problem?


texts: like finding errors overlap on the same layer in postgis?

texts:
After meeting the overlay erroneous   polygons in the same layer, I used the following query:

select gid from limite_localidad where not ST_Overlaps(the_geom,the_geom)

I have to leave 108 records error but brings me all records.
if I make the same query in qgis with topology checker, it does well giving out 108 errors, but with PostGIS I don't know how to find these errors.


texts: Clicking on data point returns "There has been an error..."

texts:
My CartoDB map has different data shown for hovering and clicking. It works fine when I'm logged in, but other folks can't see the clicking data. CartoDB thinks for a few seconds, then returns "There has been an error..."

I'm doing this map totally with the GUI, not programming anything directly. The map has two layers, one of county population and the other of about 5000 schools, and it's the latter that's the problem.


texts:
I’ve using Mapbox Studio to create a web map that’ll allow a user to interact with some basic parcel data.  I’m looking to add a basic spatial query tool for the user so that they can apply spatial filters such as: Parcel_Size >=’5000’ or Zone = ‘Residential’. Here’s an example of ArcExplorers Query Tool which has the functionality I'm looking for.

I'm new at web mapping and assuming I'm going to need to incorporate some JavaScript in order to make this work.

Does such tool exist that will allow me to add it to my Mapbox map?

I checked out Mapbox JS, but the examples I'm seeing are a little to basic (toggle buttons, search boxes). Of course I could be completely missing something here.


texts: `gdal_calc.py` only outputs as Int16 even though set to Float32

texts:
The following code only outputs integer even though set to float32.

'gdal_calc.py -A input_raster.tif --A_band 1 --outfile output_raster.tif --calc "0.005 * A.astype(float)" --NoDataValue -32767 --type Float32'

The input_raster.tif is Float32.

I have followed all the solutions here but it is not working.


texts: QGIS - Converting GEOPDF to MBTILE gives error "band 1 has no colour table"

texts:
I'm trying to convert GEOPDF files to MBTiles using QGIS on mac. I'm using GDAL version 1.11.4, QGIS version 2.8.9. I am doing this by raster -> translate -> input layer GEOPDF file, output file MBTile file. I have tried setting expand to Grey, RGB and RGBA but they all give the same error. "band 1 has no colour table".

I have installed the GEOPDF Plugin for GDAL.

I have found questions with similar issues and they all point to using "1-Band Raster Colour Table" Plugin but this does not work with the current version of QGIS.

I tried reverting back to earlier version of QGIS but then I couldn't get QGIS to read GEOPDF files.

I can use QGIS to export geotiff files from the PDF files and then use Map tiler to get the MBTile but it would be much easier if I could just do it all on QGIS.


texts: Spectral signature and its standard deviation using R

texts:
I am trying to plot a spectral signature using R similar to the one shown in the image. This plot was obtained using Semiautomatic Classification Plugin-QGIS. Half-transparent regions show the range of values, however I would like to represent the standard deviation.

I extracted the pixel values in each band and its class (land cover). My current data structure is a DataFrame like the following:
red <- c(152, 28, 34, 36, 139, 144, 29, 153)
green <- c(23, 78, 69, 68, 22, 18, 76, 26)
blue <- c(139, 18, 14, 29, 137, 162, 19, 149)
class <- c(1, 2, 2, 2, 1, 1, 2, 1)

df <- data.frame(blue, green, red, class)
df$class <- as.factor(df$class)
names(df) <- c("430nm", "530nm", "620nm", "class")

# my data structure
df

> df
    430   530   620 class
1   139    23   152     1
2    18    78    28     2
3    14    69    34     2
4    29    68    36     2
5   137    22   139     1
6   162    18   144     1
7    19    76    29     2
8   149    26   153     1

I would like to show the standard deviation instead of the range and keep the scale of the X axis and the dashed line based on the wavelength of each band.


texts: Creating DEM from contours in GRASS QGIS?

texts:
GRASS GIS Plugin - there are several modules you can use (v.surf.* and r.surf.*). You need convert your shapefile into a GRASS database to use the GRASS modules in QGIS.

I tried this option but its not processing.
what should i do?
Is their any other way in SAGA QGIS or something else please let me know.


texts: Creating Minimum Spanning Tree based on street network using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.0 SP2 (ArcView License) and Network Analyst.

Background:
I try to create a network of fiber optic cables to specific points along a street network.

Here is what I want to do:
I want to create a Minimum Spanning Tree based on a street network to give me the shortest path between a couple of fiber optic access points. The MST should start at a fiber optic distribution point and "span out" to the access points.

What I tried so far:
I tried to manage it with the Network Analyst and the 'Closest Facility' analysis. The problem is that it handles each connection between a distribution point and a access point as a independent solution (end up in a single route) - so its not the MST.
I am not sure if something like that is possible with a geometric network (with sinks and sources) and a flow direction.


texts: Extrude building footprints using raster (DSM) values in ArcGIS Pro?

texts: Is there a way to extrude a buildings footprint shapefile into 3D buildings using values from the area's DSM in ArcGIS Pro?


texts: Remove all the labels from my OpenMapTiles Server

texts:
I have setup my own OpenMapTiles Server. I have downloaded the whole Asia region. I'm using Leaflet to fetch the tiles.
I want to remove all labels from the map. I can make changes on the server or my JavaScript code but I don't know where to start.


texts: How to update length property of feature class field in ArcGIS Desktop?

texts:
I have a feature class with a string field that I want to update the length from 10 to 25.

Is there a way or a tool to update this property without having to create a new field and using field calculator to move over the records from old to new field?

This script below does what I want, the only thing is it pushes the updated field to the end of the table (not keeping the original field order).

http://arcscripts.esri.com/details.asp?dbid=16503

Another option that seemed promissing is to use the Feature Class to Feature Class tool (access this tool by right clicking a layer in ArcCatolog and selecting Export>to   Geodatabase single option) .  It does create a new layer, however you can update field names and property information.  It also appears that you can change the field order using the "Move Input Field Up/Down" arrows, but they seem to not work.

There has to be a solution to just edit field properties without having to create a new field.


texts: integrating spatialite and php (sqlite3) on windows platform

texts:
anyone can give me an explanation about how to integrate spatialite with php(sqlite3)-windows platforms.

I have tried to refer from this link http://www.gaia-gis.it/spatialite-2.4.0-4/spatialite-cookbook/html/php.html

But I do not understand is why the error message appears as follows.

"Warning SQLite3::loadextension(), Not supported in multithreaded Web Servers"

i use libspatialite-2.dll


texts: GDAL vsi (vsis3, vsizip, vsicurl) usage inside Geoserver as URL

texts:
I am trying to use vsi (vsis3, vsizip, vsicurl) from GDAL extension to read file directly from S3.
i have already loaded GDAL Native Library (version: 2.4.4) into Geoserver.
i can not use /vsis3/bucket/key as URL while i am trying to add new store in Geoserver.
do you know how can i use vsi (vsis3, vsizip, vsicurl) inside Geoserver as URL?


texts:
The link below provides a code snippet for editing an ArcSDE table with ArcObjects.

Are there any other or newer code examples out there of how to edit an ArcSDE table by using ArcObjects?

Edit SDE geodatabase table w/ ArcObjects


texts: How to display a popup for two points which have same longitude and latitude?

texts:
I have a vector layer with lots of POIs/Points. Some of the Points have same coordinates(Longitude,Latitude).

For Example:-

POINT A : {id:1,name:A,Geom:POINT(73.1545, 13.126)}

POINT B : {id:2,name:B,Geom:POINT(73.1545, 13.126)}

Current Situation:-

I can only show popup with any one of the above Point's details.
User can't able to see both POINTS on the map.

My Queries:-

How to display popup with entire details of above two POINTS?
How can we present/show these POINTS on the map?

This the select control I am using.

    var controls={};
    controls["select"] = new OpenLayers.Control.SelectFeature(LayerName,
                 {
                    multiple: false, hover: false,
                    onSelect: onFeatureSelect,
                    onUnselect: onFeatureUnselect}
                 );
    map.addControl(controls["select"]);
    controls["select"].activate();

On Selecting a feature below function will be called.

function onFeatureSelect(feature){
    var h = "", p=feature.attributes;
        h = "
";
        for(var q in p){
            h += "
"+q+":
" + p[q]+"

";
        }
        h   +="
";
        popup = new OpenLayers.Popup.FramedCloud("chicken",
                    feature.geometry.getBounds().getCenterLonLat(),
                    null,
                    h,
                    null, true, onPopupClose);
        popup.backgroundColor='#FFFFFF';
        popup.autoSize=true;
        popup.maxSize=new OpenLayers.Size(500,400);
        popup.setBorder('1px');
        popup.panMapIfOutOfView=true;
    feature.popup = popup;
        map.addPopup(popup);
}

function onFeatureUnselect(feature){
    if(feature.popup){
        map.removePopup(feature.popup);
        feature.popup.destroy();
        feature.popup = null;
    }
}


texts: Selection/Deletion of all smaller polygon/s in a SHP - Python solution?

texts:
ArcMap 10.2 user, with Spatial Analyst only. No ArcServer or of sort.

I have abt 500+ polygon shapefiles, each containing 2 or more polygons.
All of them has a NAME column and the rest of the columns may or may not be the same with each other, and may or may not has a AREA column.

My objective is to delete all the smaller polygons in each SHP, leaving only the largest one.

I know Definitive Query will not work on this as the largest/smaller polygon definition is rather relative. (Largest polygon in one SHP could be the smallest in another).

As I am hopeless with Python (still reading Dummies for Python), can anyone advise how do I go abt doing this in Python? Maybe I would need to create a AREA column for all SHP first with Python?

Another thing is the 500+ SHP is not the end, we are still adding to that folder on a weekly basis. So I would need something that we can run on with SHP consisting one or more polygons in the future.

Thanks for reading this.


texts: Vector layer does not appear with OpenLayers

texts:
When I import the vector it simply does not show in the map. I tried with a raster and it worked fine. I tried now with both KML and GeoJSON and the same thing happen, I can see that it loaded in the Chrome Dev Tool and I can also see the data from there, but not in the map. I using the Vector Layer and the Image Vector Layer as a base.

The projection is EPSG:31982, UTM 22 S - Sirgas 2000

I'm using OpenLayers v5.3

I tried some things that I found in examples online, but none seems to work.
Other questions related that did not helped me:

Vector layer not displaying in OpenLayers
OpenLayers does not display vector tiles in their origin projection
https://stackoverflow.com/questions/28688466/vector-layer-is-not-showing-in-openlayers-3
https://stackoverflow.com/questions/49769770/openlayers-4-wont-display-vector-layer
https://stackoverflow.com/questions/18124081/openlayers-vector-layer-not-showing-on-the-whole-screen
https://stackoverflow.com/questions/31879232/features-added-to-vector-layer-are-not-visible
https://stackoverflow.com/questions/13491164/openlayers-vector-features-layer-invisible-do-not-render
OpenLayers 4 does not display vector layer from geoJSON
Layer does not displaying
Not displaying any layers in OpenLayers
Openlayers 3 Circle Vectors Not Visible
https://stackoverflow.com/questions/49769770/openlayers-4-wont-display-vector-layer

Sample code I'm using:

proj4.defs('EPSG:31982', '+proj=utm +zone=22 +south +ellps=GRS80 +units=m +no_defs');
register(proj4);

var map = new Map({
  layers: [
    new TileLayer({
      source: new OSM()
    }),
    new VectorLayer({
      renderMode: 'image',
      source: new VectorSource({
        url: 'img/estradas.geojson',
        format: new GeoJSON(),
        projection: 'EPSG:31982'
      }),
      style: function(feature) {
        style.getText().setText(feature.get('name'));
        return style;
      }
    })
  ],
  target: 'map',
  view: new View({
    center: transform(getCenter(imageExtent), 'EPSG:31982', 'EPSG:3857'),
    zoom: 13
  })
});


texts: Positional accuracy assessment methods in qgis

texts:
I'm doing an quality assessment of openstreetmap's street data. As a comparism dataset I have to use an administrative dataset, which is highly accurate and comes in polygons. Here's a sneak peek:

So far I've examined positional accuracy via goodchild and hunter's buffer method. But that's not enough content and I'm struggling to find more evaluation methods I could build some models with. I came across hausdorff distance, but a (required) conversion from polygons to polylines didn't work out well.

Any ideas?


texts: Where can I download historical data from OpenStreetMap? I am interested in downloading an OSM .pbf file of a particular city in previous years (say, from 2008 and 2010).


texts: button background image from databse to external css

texts:
I am adding buttons dynamically based on the data in the database and i want to dispaly the buttons with images(each button have differnet images). images path is also defined in the database. my doubt is how to pass this button images from database to css.

my code is like this

var button1 = new OpenLayers.Control.Button({title:arraytypes[i], displayClass: "button1Display", trigger: selectmarker});

"button1Display" is the class in the extenal style sheet

this is my style sheet

.button1DisplayItemInactive
{

  background-image: url("../smartcity/images/marker.png");

  background-size: 50px 52px;

  background-repeat: no-repeat;

}

in the above code background-image is hardcoded but i want to change it dynamically (from database).


texts: Adding ArcGIS-FeatureServer-Layer to QGIS using PyQGIS?

texts:
I have a password protected ArcGIS-FeatureServer Service, the service provides 25 layers.

Using "Add Layer --> Add ArcGIS-FeatureServer-Layer" works like it should.

Is there any way to add these kind of layers with PyQGIS? I could not find a documentation or example anywhere.

How to Add a normal WFS layer with PyQGIS is well documented.

For ArcGIS MapServer it is super-easy:

iface.addRasterLayer("url='my_featuere_layer_url' layer='0'", "mylayer", "arcgismapserver")

But the provider "arcgismapserver" is not documented anywhere, just found an example after some search, so I expect that there is a similar solution for ArcGIS-FeatureServer.


texts: How to store a drawn feature in OpenLayers?

texts: I have found a draw tool in OpenLayers-3 which allows us to draw on a map (see http://openlayers.org/en/v3.2.0/examples/draw-features.html). Is there an API in openlayers to store the drawn feature in geojson format?


texts: Increase number of parallel caching tasks for GeoWebCache?

texts:
By default GWC has limitation of 16 parallel caching tasks.

Is it possible to increase this limitation?


texts:
Do you know how to call MRT Swath from IDL?

I would like to work with Modis L1B (MOD02HKM), read it, reproject, some processing and writing the results.


texts:
I am using CesiumJS to create point in a map using World Terrain to get the elevation of my points.

The point in Cesium are using WSG84. When I calculate the altitude using sampleTerrain I retrieve an height above the reference ellipsoid

To convert that height to MSL, I have seen on Internet Here or Here that require a geoid to compute the conversion from WGS84 point to MSL.

  The geoid height above the ellipsoid (N) is the difference between the
  ellipsoidal height and orthometric (geoid) height.

  N = He − Ho

  Where:

  N = geoid height above the ellipsoid
  He = ellipsoidal height
  Ho = orthometric (geoid) height

I am using a package that calculate a geoid using EGM2008 for a given long/lat.

So I dispose of :

A Point A in WSG84 with an height above the ellipsoid
A Point B in EGM2008 with a geoid height

I tried to do MSL = A - B but my value are always wrong.

I used those two calculator :

unavco
daftlogic

When using those two tools, when I do GPS - GEOID = MSL it does work, but the value I get from Cesium are not the same (Isn't GPS using WSG84 ? )

 GPS  |  CESIUM       | GEOID      |  DIFF (GPS / CESIUM)

669   | 705           | 49.421     | +36
52    | 88            | 37.162     | +36
51    | 36            | -15.992    | -15
0     | -86           | -98        | +86

So my question is : Using Cesium World Terrain and apparently WSG84 altitude above an ellipsoid, how do I retrieve the altitude of the point in MSL?


texts: Value of Raster when we combine 2 raster. (New value not = with old value)

texts: I have two raster in ArcGis (Number 1 is water level, Number 2 is water depth). I want to combine it to only one raster with two field value. I used Combine Tool in toolbox in ArcGis. After that I got one raster but two field contain value of raster Number 1 and raster Number 2 not = with value of Number 1 and Number 2. Please help me know why and which what order way to have good result.


texts: gdal_rasterize and keep/create the attribute table

texts:
I want to rasterize my polygon by a given field, and keep/create the attribute table.
Now, as a raster, the attribute table will contain a row by each unique value (category) in my raster, and no an entry by each polygon, as vector's attribute tables. ArcGIS rasterization creates an *.var.dbf (or so) table, counting the pixels per category. It is what Im looking to do.
Any idea of how do this using preferebly gdal_rasterize?


texts: "QGIS-BIN.EXE has stopped working" error after QGIS run

texts:
After installing QGIS 2.18.7.1 on a 32-bit Windows 7 machine, when I run QGIS I have the following message: "QGIS-BIN.EXE has stopped working"

Do you have the solution to the problem?

I tested version 2.18.4, 2.18.6 ==> IDEM same problem.


texts: Mann Kendall test on snow cover images in Google Earth Engine

texts:
I am trying perform Mann Kendall test in Google Earth Engine, but output show zero image. I am attaching the code link. https://code.earthengine.google.com/a6dd42809be0452450990cec70de71b0

dataset = ee.ImageCollection('MODIS/006/MOD09A1')
              .filter(ee.Filter.date('2001-01-01', '2002-12-31'))
              .map(function(img) {return img.clip(table)});

var cloud_free = dataset.map(
  function(img) {
  var qa = img.select('StateQA')
  var cloudState = bitwiseExtract(qa, 0, 1)
  var cloudShadowState = bitwiseExtract(qa, 2)
  var cirrusState = bitwiseExtract(qa, 8, 9)
  var mask = cloudState.eq(0) // Clear
  .and(cloudShadowState.eq(0)) // No cloud shadow
  .and(cirrusState.eq(0)) // No cirrus

  var maskedImage = img.updateMask(mask);
  return maskedImage.addBands(mask.rename('mask'));
  })

 function bitwiseExtract(value, fromBit, toBit) {
  if (toBit === undefined) toBit = fromBit
  var maskSize = ee.Number(1).add(toBit).subtract(fromBit)
 var mask = ee.Number(1).leftShift(maskSize).subtract(1)
 return value.rightShift(fromBit).bitwiseAnd(mask)
 }

 var snow_cover =              cloud_free.map(function(img1)
              {
                 var ndsi = img1.normalizedDifference(['sur_refl_b04','sur_refl_b06']);
                var evi1 = ndsi.expression('(thr>0.4) ? 1 : 0',
                                           {  'thr' : ndsi  });
                 return evi1.clip(table).updateMask(img1.select('mask'));
              });

   var afterFilter = ee.Filter.lessThan({
   leftField: 'system:time_start',
   rightField: 'system:time_start'
   });

    var joined = ee.ImageCollection(ee.Join.saveAll('after').apply({
    primary: snow_cover,
    secondary: snow_cover,
    condition: afterFilter
     }));

     var sign = function(i, j) { // i and j are images
     return ee.Image(j).neq(i) // Zero case
     .multiply(ee.Image(j).subtract(i).clamp(-1, 1)).int();
     };

     var kendall = ee.ImageCollection(joined.map(function(current) {
     var afterCollection = ee.ImageCollection.fromImages(current.get('after'));
       return afterCollection.map(function(image) {
        // The unmask is to prevent accumulation of masked pixels that
      // result from the undefined case of when either current or image
      // is masked.  It won't affect the sum, since it's unmasked to zero.
       return ee.Image(sign(current, image)).unmask(0);
      });
     // Set parallelScale to avoid User memory limit exceeded.
       }).flatten()).reduce('sum', 2);

      var palette = ['red', 'white', 'green'];
    // Stretch this as necessary.
     Map.addLayer(kendall, {palette: palette}, 'kendall');
     print(joined)


texts: Using field Calculator for time in ArcGIS 10.2

texts: I have a  new field in a shapefile with the property "Date". I want to fill it with "9:00:00 AM". I tried several ways of using field calculator, but it gives date instead of time. How could I do that?


texts:
I am using ArcGIS 10.2, and previously using Oracle spatial as back end.

Now the data providers have shifted to SQL providing core data using FTP which is being staged at SQL server.

We need to import the spatial data from SQl server on daily basis and place in Oracle so it doesn't effect service delivery.

Which would be the best way to add SQL Server connection through ArcGIS desktop?


texts: How long is a GPS epoch in microseconds?

texts:
My GPS is reporting dates from the previous epoch.  How many microseconds should I add to the date to get the correct value?

Since the 2019 GPS Week Number Rollover, my old GPS logger (QStarz BT-Q1000X) has started reporting that the current GPS date is in 1999.  I think I can fix it by re-writing the logging software to add a constant value to the timestamps.  But what constant?  The software is written in Java, which measures time as the number of microseconds since midnight at the start of 1 Jan 1970, so I need a value in microseconds.

From the answer to What is a GPS epoch? I understand that GPS time is measured in weeks and seconds from 6 Jan 1980.  In this sense, the length of the GPS epoch is the time taken for the week number to roll over.  This is 1024 weeks (as it is stored as a 10-bit number) or about 19.7 years.  Hence the second GPS epoch started on 22 August 1999, and the third on 7 April 2019.

The simple calculation is that I should add on 1024 weeks * 7 days * 24 hours * 60 minutes * 60 seconds * 1000 = 619 315 200 000 ms.  But is this right?  Do I need to account for leap seconds, non-integer numbers of weeks/year or hours/day, etc?  In short, is there anything else non-obvious I need to take into account?

For my application, I don't need accuracy better than to the nearest second, and I realise that there is a difference between GPS time and UTC.  What I'm after is values that are reasonably accurate with respect to current GPS time.

More relevant background:

Reference for Java's Date class in Android.
Ticks between UNIX epoch and GPS epoch (discusses the offset between the start of UNIX time and GPS time).


texts: OpenLayers 4 popup from multiple layers with different source

texts:
When I click on the map I need to create a popup and query the 3 layers (at the time I'm writing I have just 2 but they will be 3).

Here my code:

  amfeatlayer = new ol.layer.Tile({
    extent: ext,
    preload: Infinity,
    source: new ol.source.TileWMS({
      url: amfixed.wmsfeat,
      params: {
        'LAYERS': rtctwlayers,
        ...,
        format: 'image/png',
        'TILED': true
      },
      ratio: 1,
      gutter: 40,
      isBaseLayer: false,
      serverType: 'mapserver',
      crossOrigin: 'anonymous',
      projection: 'EPSG:27700',
      tileGrid: new ol.tilegrid.TileGrid({
        resolutions: resolutions,
        extent: ext,
        tileSize: [128, 128]
      })
    })
  });

  pmfeatlayer = new ol.layer.Tile({
    extent: ext,
    preload: Infinity,
    source: new ol.source.TileWMS({
      url: pmfixed.wmsfeat,
      params: {
        'LAYERS': sQLayers,
        MYORDERS: twPmSettings.defsel,
        DTO: sDTo,
        format: 'image/png'
      },
      ratio: 1,
      gutter: 40,
      serverType: 'mapserver',
      crossOrigin: 'anonymous',
      projection: 'EPSG:27700',
      tileGrid: new ol.tilegrid.TileGrid({
        resolutions: resolutions,
        extent: ext,
        tileSize: [128, 128]
      })
    })
  });

and here the code I use to get a popup:

map.on('singleclick', function(evt) {

var content = document.getElementById('popup-content');
content.innerHTML = '';
var viewResolution = map.getView().getResolution();
var coordinate = evt.coordinate;

url = pmfeatlayer.getSource().getGetFeatureInfoUrl(coordinate,
viewResolution, projection, {
  'INFO_FORMAT': 'text/html'
});

if (url) {
  var xhttp = new XMLHttpRequest();
  xhttp.open("GET", url, true);
  xhttp.send();
  xhttp.onreadystatechange = function(aEvt) {
    if (xhttp.readyState === 4 && xhttp.status === 200) {
....

How can I query all the available layers?
As you can see I have different source for every layer, I tried to combine the 2 URLs but it doesn't work properly (just to check if I can have the URL of both and then split "---" and pass each URL in the function)

URL = pmfeatlayer...
URL + "---" + amfeatlayer..

but it still printing out the first URL


texts: Problem with writing the srsName into gml output in OpenLayers

texts:
I am trying to write out my geometry in gml format after it is drawn in OpenLayers. Actually it is working fine, I get the gml like this:

<gml:Polygon>
<gml:outerBoundaryIs>
<gml:LinearRing>
<gml:coordinates decimal="." cs="," ts=" ">8.62975388121127,46.34818691143132 8.761589818711352,46.590322794368646 9.124138646836665,46.43911409951752 8.805535131211348,46.24951148519901 8.62975388121127,46.34818691143132</gml:coordinates>
</gml:LinearRing>
</gml:outerBoundaryIs>
</gml:Polygon>

But unfortunately I need it to have the coordinate system (srsName) defined in it, like this:

<gml:Polygon srsName="EPSG:4326" >
<gml:outerBoundaryIs>
<gml:LinearRing>
<gml:coordinates decimal="." cs="," ts=" ">8.62975388121127,46.34818691143132 8.761589818711352,46.590322794368646 9.124138646836665,46.43911409951752 8.805535131211348,46.24951148519901 8.62975388121127,46.34818691143132</gml:coordinates>
</gml:LinearRing>
</gml:outerBoundaryIs>
</gml:Polygon>

Here is my code:

draw.events.register('featureadded', this, function(evt) {
   var feat1 = vectors.features;
   var gml = new OpenLayers.Format.GML({
                  format: OpenLayers.Format.GML.v3,
                  formatOptions: {
                       srsName: "urn:x-ogc:def:crs:EPSG:4326",
                  }
             });
   for (var i=0; i<feat1.length; i++) {
   var geomtransformed = feat1[i].geometry.transform(new OpenLayers.Projection("EPSG:3857"), new OpenLayers.Projection("EPSG:4326"));
   geomout = gml.buildGeometry.polygon.apply(gml,[geomtransformed]);
   console.log(geomout)
   }
});

This is working but the srsName does not appear in the gml.

So could you please advise me if it is possible to transform vector feature to gml with writing the srsName into the gml, and if yes, how?


texts: NMEA to CSV with timestamp using gpsbabel


texts: QGIS Valmiera, Labels only show for part of polygons

texts:
Using 2.2.0-Valmiera on Ubuntu 13.10

I have for Croatia 21 Counties as polygons.
Labels show only for 12 counties: here I see that no labels are shown for multipart features.
So I checked "Show labels for multipart features": No change
I unchecked "Discourage labels from covering features: No change

I tried a layer with all Municipalities cca 500 polygons.
Labels show for most polygons but not all. All multipart polygons failed to show labels.

Then I found out that this behaviour of is dependent of the scaling. At 1:25000 all labels show, also on multipart but at 1:50.000 part of the labels are not shown.
I controlled:
- Scale based visibility : unchecked so okay
- Topology: No errors
- Geometry validity: no errors

Has anyone seen this behaviour and perhaps know a solution

Thanks


texts: Using Python Lists in Calculate Field Code Block?

texts:
I'm looking for some help in explaining what exactly is occurring with the following code:

MaterialLife = "LifeMaterial(!"+Material+"!)"
codeblock2 = """def LifeMaterial(mat):
    material = {"ABS":90,"AC":100,"BRICK":100,"CI":70,"CLAY":100,
          "CONC":100,"CONPP":100,
          "CONR":100,"CONX":100,
          "CORI":70,"CORM":90,
          "CSP":70,"CSPA":90, "CSPH":70,
          "CSPP":70,"DI":60,"PB":80,
          "PE":80,"PP":80,"PVC":100,
          "PVCC":100,"ST":80,"STC":100,
          "STONE":100,"TP":40,"PFR":90,
          " ":90,"UNK":90,None:90}

    life = material [mat]
    return life"""

arcpy.CalculateField_management(ins_pipe, MaterialLifetime, MaterialLife, "PYTHON_9.3", codeblock2)

I know the basic structure in that there's an expression (MaterialLife) and a codeblock (codeblock2) that are used to perform a calculate field function. The input field to the expression (Material or mat in the codeblock) contains all of the different types of materials that are in the codeblock ("CONC", "CONPP" etc.). The numbers beside each of the materials are ages that are to be assigned to each material. The end result of the function is that all of the ages associated to the material types are populated in the field "MaterialLifetime" based on its corresponding Material type.

My main questions are:

How does the program know to properly go through every entry in the table even though they are not listed in the order that is laid out in the codeblock without any if statements?
How does the program know to populate the "MaterialLifetime" field with ONLY the age and not the material type as well?


texts: Why rasters of the same projections don't overlap?

texts:
I tried to search for existing questions to see if I could find answers, but I didn't.  I have two rasters as shown.  The raster with gray borders is not titling the right way.  It should have tilted and overlapped with the raster with red border.

Both rasters were projected the same way, figure attached.

What caused the non-overlap of the two rasters?


texts: GeoServer serves a TIFF file in an Azure Blob

texts: Can GeoServer publish a TIFF file stored in an Azure Blob. i.e the TIFF file is not on the same machine of GeoServer?


texts:
I've inherited some old CARIS GIS data (NTX format) and am needing to convert it into something usable while ensuring that all attributes remain intact. I do know of the NTXShape application (http://ntxshape.sourceforge.net/) but have gotten little success with that, primarily because a working copy of ArcView 3.x is rare as hen's teeth and Windows 7/8 doesn't seem to like the command line application.

Are there any other applications (preferably opensource) that may be able to undertake the conversion?


texts: The "Azimuth and Distance" plugin does not appear to be maintained any longer, and has not been transitioned to the new QGIS 2.0 plugin repository. What are the current options for entering COGO in QGIS? Does the roadmap include these tools, or is this something that will still depend on a plugin?


texts: Fixing advanced labeling for a layer in ArcGIS Pro 2.6

texts:
I opened a NetCDF file in ArcGIS Pro using Multidimensional raster layer under add data. Then I assigned some symbology to the layer using the advanced symbology option within the symbology pane. Within the ArcGIS Pro environment, the symbology was assigned fine and looks like this:

But when I publish this layer as a web service and open it in ArcGIS Online, then the symbology of the layer is displayed in this way:

How can I remove the ADVANCEDLABELS text from the symbology?


texts: Converting from unknown coordinate system into common one using fiducial points?

texts:
I have a thousand of points that I want to project in a different coordinate system. These points were generated using an old and odd coordinate system (let's call it CS) so it's not possible to use some library to convert them. Luckily for each 50 points I have 1 point in both coordinate systems (the old and odd CS system and the WGS84). So how can I convert the remaining 49 points to WGS84. I think it is matter of interpolation but I don't know the math to use.

I have the WGS84 and the CS coordinates of the red points,
I need the WGS84 coordinates of the black points.

I know how to do this in QGIS using a custom CRS string with 1 point but here I need at least 3 points because 1 point will only move the layer. I also need to scale or deform.

Is it possible to have a custom CRS string based on 3 points?

I am also interested by the math to use (I will eventually do the math in Python).


texts:
I'm tryng to use the tool gdalwarp, but in some cases there is an alert that say:

  Warning 1: Self-intersection at or near point 2647000 4441118.8085795399
  ERROR 1: Cutline polygon is invalid.

Can you tell me how to fix the input shape and solve this problem of "self-intersection"?

shape_filtered=outShapefile
warp_opts=gdal.WarpOptions(format="GTiff", cutlineDSName=shape_filtered, cropToCutline=True, dstAlpha = True )

res=gdal.Warp(output_dir+input_raster[-10:-4]+"_"+cod+".tif", input_raster, options=warp_opts)
print (res)


texts: Satellite imagery alternatives for land use/cover classifcation

texts: Apart from Landsat imagery, which other satellite imagery is good for land use / cover studies and can be obtained free of charge?


texts: Generating line segments between all points using QGIS

texts:
How do I create a vector layer with a network of line segments drawn between every point in another vector layer?
Basically, I am trying to create something similar to a hub/spoke network, except that instead of having a single hub and the other points are destinations, I would like every point to be treated as a hub with a direct connection to every other point.
I have managed to get the hub lines from one point to every other using the MMQGIS plugin "Hub Lines" module, but it seems that I would need to re-run it for each individual point.
Is there a way to just do it all in one go in QGIS?
It may be worth noting that I don't need to know the distances between the points; I'm looking for the azimuths.


texts: How to upload polygon data for appending to a geodatabase via a JS API web application

texts:
After a many hours of research, I am having a hard time determining the feasibility of an application I am trying to build.

Basic proposed functionality is: user logs in, uploads polygon (either shapefile or feature class), and is then able to attribute the feature using the new server 10.1 Feature Access capability.

The file geodatabase to be appended is on my own ArcGIS Server instance, and is connected to an Oracle enterprise DB.

I am aware that 10.1 allows a user to build a polygon right in the application window using feature access and editing, but the nature of the geometry being submitted will not allow for that. The data will have to be created in ArcMap by the app users, and then uploaded by them to a gdb on a server my organization maintains.

While I did indeed read the accepted answer to this question, I see no mention there of making edits to a geodatabase which is remote from the application user, and so I am unsure that it actually addresses my issue. I would gladly be corrected on that.

I know that at least a portion of my proposed functionality will be accomplished using a geoprocessing service (the appending of the new data to the existing gdb), but its the uploading of the data where I am stumped.

I am nearly convinced that this is not quite possible using tools available in the ArcGIS API. Does anyone think it is, and have any tips on where to look?


texts:
Based on this: Apply an SLD saved on Geoserver to OpenLayers 3 I try to set a SLD in the params of a TileWMS. The SLD is saved on Geoserver, but everything i tried failed.
The strange thing about it is, if I define the style in a normal GetMap request it works.

My main goal is that the user can decide in which style the layer should be styled with something like that:

WMS_source.updateParams({'STYLES':'name_of_style'});


texts: Error: Computation timed out, EE memory exceeded in reduce region.flatten()

texts:
I am trying to run a fairly in-depth script that reduces various environmental covariates to HUC12-level polygons and exports the information to a csv file. My script gets an Earth Engine capacity exceeded error at the following x marks.

Here is a brief outline of my script:
1.) Select variables of interest from various EE assets (n in total) x
--
2.) Map over a list of years (14 in total):
2a.) Grab annual and seasonal (xx) summary statistics
2b.) Add each summary statistic as a separate band to the year of interest
2c.) End with a list of length 14, with n bands
--
3.) Reduce the banded image over polygons of interest (in this case HUC 12 units) xxx
3a.) End up with a list of Feature Collections of length 14
3b.) Each FeatureCollection contains every polygon of interest which has been reduced
3c.) Embed client-side information into banded image via an ee.Join

The X's indicate places where an "Earth Engine Capacity Exceeded" error has occurred.
Sequentially:
x I have been able to run the analysis such that I end up with n = 28 covariates, but I am interested in > 40
xx I am interested in computing, for example, total seasonal precipitation for each year, which adds about 4x as many bands
xxx I am able to acquire a list of Images (end of step 2) with all covariates of interest, but am unable to then reduce them to my polygons of interest

A brief look at my profiler reveals that the process with the highest memory only is using about 17 megabytes.
I have tried running the script without additional print() requests, as detailed here, and have referred to this question as well, regarding the order in which I map to ImageCollections.
This link seems to suggests that the Earth Engine memory capacity exceeded may be a bug, but I am well aware that this is a computationally expensive project.
Here is a link to my script.
Edit: I have tried submitting the server request from a different network and the analysis still failed. I am wondering if there is a more efficient way to get the information that I seek
Here is a link to the EE Dev Thread


texts: PostGIS crosstab on multiple point tables

texts:
I have several point tabels, lets say 3. Tables are alike with two columns: a integer value and a geometry like:

table sand1 (value int, geom geometry)
table sand2 (value int, geom geometry)
table sand3 (value int, geom geometry)

I need to join the three tables to a new table with four columns like:
valueTab1, valueTab2, valueTab3, geometry.

The join has to be on the point geometry. The value field cannot be used since it is a measurement value. Some of the points are positioned directly on top of each other. The ST_Equals points then results in one record.

So the picture above with three layers sand1 (three records), sand2 (two records) and sand3 (three records) should produce a new table with four columns and four records.

Columns:

sand1Value, sand2Value, sand3Value, geom

Records:

22, null, null, null, geom
null, 31, 37, null, geom
1111, null, 66, geom
11, 44, 55, geom

I have tried joining with ST_Equals(sand1.geom, sand2.geom), but ST_Equals does only take two parameters. Have tried the tablefunc extension with crosstab and multiple joins on like WHERE ST_AsText(sand1.geom) IN (ST_AsText(sand2.geom), ST_AsText(sand3.geom))


texts: Why is GeoTIFF squashed in WGS84 data frame of ArcMap?

texts:
I have a government-sourced geotiff which I have imported into ArcMap 10.4 as a raster basemap. In a TIFF viewer it displays correctly but in ArcMap it appears vertically squashed by roughly 20%.

The data frame coordinate system is GCS_WGS_1984 WKID: 4326. And here is some data provided to me along with the tiff file:

Panel Limits:      Lat  -34 5 55.812    to  -33 23 59.347
                       Lon  151 7 58.170    to  151 39 42.778

Referencing:       Geographic Lat/Long
Spheroid:          WGS84
Horizontal Datum:  WGS 84
Scale:             1:75000

Image Size:        7258 cols x 10867 rows
Image Resolution:  254 dots-per-inch
Number of Colours: 27

TIFF Version #:    42
GeoTIFF Version #: 1
TIFF Compression:  Packbits

Why does it appear squashed, and how can I fix it?

I found this but it did not seem to help: https://geonet.esri.com/thread/34535. Unable to independently verify the geotiff is WGS84. I am running Windows 10 Pro under Parallels on a MacBook Pro.

What is stranger, I also have a KML map of the same area and when I alternate between the two maps in ArcMap, by disabling one of them, I can see that they are exactly coincident. In other words, both are squashed, by the same amount.

The only thought I had was that one degree of longitude, in metres, varies with the cosine of the latitude, and coincidentally 1 - cos (33.5) is about 20%. Could that have something to do with it? Clutching at straws here.


texts: Using QgsExpression and expression strings for postal address search

texts:
I think this question is more about using QGIS expression strings. I have a large shapefile of postal addresses in an area. And I am trying to write a script where the user inputs an address and the script returns the feature using QgsFeatureRequest. In the shapefile is an address field BEZ which contains the concatenated full address. This is what I have so far

#address to lookup
address_string = address_string.upper()
address_split = address_string.split(' ')
expr = ""

#expression string
for i, add in enumerate(address_split):
    expr += "\"BEZ\" LIKE \'%{}%\'".format(add)
    if i<len(address_split)-1:
        expr += " AND "

address_ftr = addresses_lyr.getFeatures(QgsFeatureRequest(QgsExpression(expr)))

The expression string works for street name and type, but it doesn't exactly work for the street number. Here are some cases where it doesn't work:

Looking up "42 Main St" returns "42 Main St" and "142 Main St"
There are cases where the address is an amalgamated lot so the street number goes like this "142-143 Main St", which would also be returned when you search "42 Main St"


texts: Is there a method to create watershed on google earth?

texts:
I want to create a watershed boundary for a lake only using google earth without using other source of data like srtm etc.

Since we get the height of the location can we use it somehow and create a watershed boundary with effort than digitizing?


texts:
I am fairly new to the GIS and wanted to ask something which can be very preliminary or even worse totally invalid.

How do I create a geom randomly centered around or close to a given point anywhere.

i.e. for example

{"type":"Feature","id":"01","properties":
{"name":"Alabama","density":94.65},"geometry":{"type":"Polygon",
"coordinates":[[[-87.359296,35.00118], [-85.606675,34.984749],
.....,[-88.098683,34.891641],[-88.202745,34.995703],[-87.359296,35.00118]]]}},

so the geometry field value in the above geojson is for Alabama region.

If I were to create random complex polygon geom (with at least 20 different points) on map through code, how would I do it?

Basically I want to generate a random polygon around a point and this polygon can't be simple shape (i.e. quadrilateral or so)


texts: Algorithm "qgis:difference" doesn't work in pyQGIS 3

texts:
I try to use native algorithm difference of QGIS. But it seams it doesn't work like the others :

parameters = {'INPUT': polygons['OUTPUT'], #from temporary layer created before
              'OVERLAY': layer, #read in the beginning with QgsVectorLayer
              'OUTPUT': 'memory:'}
differenced = processing.run('qgis:difference', parameters)

I get the following error :

_core.QgsProcessingException: There were errors executing the algorithm.

Both input and overlay layers have the same geometry. Any idea ?


texts:
Question

How do I convert an SVG file to GeoJSON format?

Desired Outcome:

I want to make a choropleth map using something like D3.js, or any of the other libraries out there.

Challenge:

I have an SVG file that I was able to obtain from public domain of the administrative regions for a particular province of a particular country.  I believe I need to convert that to GeoJSON before being able to use it in the libraries that I am aware of.


texts:
I would like to create two presets in QGIS 2.18 that are displaying the same layers, but the first preset using one labels/colors and the second preset using another labels/colors.

As a workaround, I "duplicate" the layer and display the original layer for preset A and the "duplicate" for preset B, with two different labels & style.
However this can create a lot of duplicate layers in big projets.

Is there a function "preset" to define the rule-based labelling and the rule-based style?


texts: How to remove points which are too close to each other?

texts:
I have a bunch of data points, however I think that some might be too close to each other due to sampling error and could therefore have an impact on my subsequent analyses.

How can I remove data points which are less than 200 meters from each other?

I'm using ArcGIS 10.2


texts:
I have bounding box coordinates of my shapefile
(fetched from layer.GetExtent() Python with GDAL/OGR):

(-38.6102467739459, -38.017601026186576, 33.01563382506067, 33.624945228091406)

So I assume it's SW and NE points of the rectangle.

I would like to find out the center of this rectangle (in lat/lon).
How should I proceed?


texts: What page size to insert a map into a report?

texts:
50% of my job is creating map figures in ArcGIS, most of which end up being inserted in Word documents. We use A4 MXD templates, but the maps end up much smaller when inserted into a report. So it seems a bit silly using A4 templates, when the figures end up being resized.
Taking into account 'standard' margins in word processing software (and using A4, not ANSI paper size), what is the optimal page size in which to export a map? They may also have a small caption beneath. I may have just answered my own question, after I go and measure standard margin sizes... But hoping there may be some useful insights from fellow industry colleagues.


texts: Cannot get geometry column name from geopandas

texts: I tried calling gdf.geometry.name and I always get back "geometry" as the geometry column name. Using GDAL's ogrinfo, I get the actual name which is "geom". Is there another way of getting the correct geometry column name without parsing the output from ogr?


texts: how to convert shp and gis line and raster data to autocad (dwg or dxf) in QGIS or GRASS

texts:
I am trying to find out which is best free GIS program with ability to convert shp and line and raster files to files suitable for autocad (e.g. dwg and dxf).

I have QGIS and GRASS.

I get an error message that QGIS cannot save as dxf although it is a format in the options.


texts: Translating vector by raster DEM data in QGIS?

texts:
I'm trying to follow the set of instructions here: https://twitter.com/imhere_asia/status/1023055838787690496

[EDIT]
Pasting the instructions here for reference so people can avoid twitter.

  Add a DEM raster into your project
  Open the processing panel, and select the "Create grid" algorithm
  You'll want to create a point grid, with narrow horizontal spacing (roughly matching the cell size of you DEM raster) and larger
  vertical spacing)
  Select the "Translate" algorithm
  Use the above-created point grid as your input layer
  The magic happens here: set the y-axis offset distance value to be the following data-defined expression:
  "raster_value('dem_layer',1,$geometry)
  Run the algorithm
  Lastly, select the "Points to path" algorithm
  Use the above-created translated layer as input layer
  Use the "left" field as order field and "top" field as group field
  Run the algorithm, et voila.

The end goal of this is to end up with a map corresponding to the DEM that looks sort of like the classic joy division cover (see image above from twitter). The spot where I'm getting stuck is step 6.

[/EDIT]

However, when I translate the grid by raster_value('my_dem',1,$geometry), the algorithm runs, but the new layer is identical. The dem layer is an SRTM GeoTIFF. What am I doing wrong?


texts: Duplicated place names when rendering tiles with mapnik and mod_tile

texts: So I am trying to render all the tiles for the state of Florida with the file generate_tiles.py that came with mapnik. The tiles are rendering, but erroneously including multiple city names (see below). Does anyone know how to fix this?


texts: List of Offensive Symbols/ Colours for Cartography

texts:
I was at a data visualization workshop a few days ago and someone there asked if there are any symbols/colours that might offend people when making charts, the presenter didn't really understand the answer so I helped to explain how in cartography using certain colours to represent certain groups could offend.

For example using yellow to represent East Asian's, Red to represent Native Americans, or Red to represent basically any minority (since red has a connotation of being negative/bad/dangerous)

The only symbols I can think of that have negative connotations are yellow six pointed stars to identify Jewish people, and pink triangles to identify gay people as these were used in Nazi Germany.

I realized I only ever learned this from working with other researchers and was wondering if there is some comprehensive list out there for these types of colours/symbols.


texts: Is it possible to fix a vector layer after I badly georeferenced a map?

texts:
Sorry in advance if I missed another question on this topic. I searched for everything I could think of, but I didn't find anything with a solution to my problem.

I am quite new to GIS tools and recently georeferenced a map of Germany. I didn't notice that I did a poor job on it though and although in some parts it is correct other's are way off. The problem I have now is that I created a vector layer with line objects based on the wrong georeferencing I did on the map.

I know I can fix the original georeferencing quite easily, but is there a way though to fix the vector layer more or less automatically? Is there any way I can somehow adjust the position of the objects in my vector layer according to the new georeferencing of that map? I would very much like to avoid having to manually create all those line objects again, but don't know how I could possibly fix it otherwise.

Thanks a lot to anyone who can help me with this, it's much appreciated! :)

Edit: I forgot to mention, I am using Quantum GIS 1.8, but if it was possible with any other Freeware I'd be open to try that


texts: Clipping by looping through list using ArcPy

texts:
I'm having trouble creating a looping function that looks through a geodatabase, clips each feature class to a "study_area" polygon, and puts the clipped files in a separate geodatabase. I keep getting "Could not clip feature classes". The features can all be clipped individually but trying to run this loop doesn't work.
What am I missing here?
I tried to search for similar questions but only found answers regarding raster clips. Code below.
from arcpy import env

# setting my source and target geodatabase
SourceWorkspace = r"E:\UFD Commercial Corridor\Roosevelt.gdb"
TargetWorkspace = r"E:\UFD Commercial Corridor\Clipped.gdb"

#setting my clip feature (which is a polygon)
clipfeature = r"E:\UFD Commercial Corridor\Clip area.gdb\study_area"

env.workspace = SourceWorkspace

try:
    fclist = arcpy.ListFeatureClasses()

    for fc in fcList:
        outClipFeatureClass = targetWorkspace + "//Studyarea_" + fc

        arcpy.Clip_analysis(fc,clipfeature,outClipFeatureClass)

except:
    arcpy.AddError("Could not clip feature classes")
    print ("could not clip feature classes")
    print arcpy.GetMessages()


texts: Version editing for Postgresql 9.2.4 / postGIS 2

texts:
Background:
I am "building" a GIS with PostgreSQL 9.2.4 / PostGIS 2 as the database and QGIS 1.8.0 as the geographic application. There will be multiple users accessing, editing and viewing the data at the same time.

Question:

Is version editing automatically enforced in PostgreSQL 9.2.4 / PostGIS 2?
Is it possible to enable/enforce version editing and if yes how?

I have had a look at pgVersion for PostgreSQL 9.2.4 / PostGIS 2 and the QGIS pgVersion plugin. However, if I understand correctly (reading the help files and testing it), pgVersion only allows editing of cell values of existing columns and also no new geographic features can be added.


texts: Let QGIS 3.0 Processing algorithm output a VectorLayer loaded via the delimited text provider

texts:
Following up on this question, I need to be able to setup the QgsVectorLayer loaded via the delimitedtext provider as my algorithm output.

Note that I am still working from this template.
For lack of a better idea, I have defined (in initAlgorithm) the output like so:

self.addParameter(
    QgsProcessingParameterVectorDestination(
        self.OUTPUT,
        self.tr('Output'),
        type = QgsProcessing.TypeVectorAnyGeometry
    )
)

And my Algorithm currently looks like this:

def processAlgorithm(self, parameters, context, feedback):

    ## Define some parameters and functions to handle the loaded file.

    # Getting parameters values to local variables
    crs = self.parameterAsCrs(parameters, self.CRS, context)
    inputFile = self.parameterAsFile(parameters, self.FILEIN, context)
    ver = self.parameterAsEnum(parameters, self.VERSION, context)

    dest_id = self.parameterAsOutputLayer(parameters, self.OUPUT, context)

    uri = ('file:///{}?'
              'type=regexp'
              '&delimiter={}'
              '&skipLines={}'
              '&useHeader=No'
              '&trimFields=Yes'
              '&xField={}'
              '&yField={}'
              '&crs={}'
              '&spatialIndex=Yes'
              '&subsetIndex=no'
              '&watchFile=Yes').format(inputFile,
                                       regex,
                                       skips,
                                       'Easting',
                                       'Northing',
                                       crs.geographicCrsAuthId())

    vl = QgsVectorLayer(uri, 'MY_LAYER',  'delimitedtext')

    # Setup Fields Aliases
    for idx, alias in enumerate(myFields):
        try:
            vl.setFieldAlias(idx, alias)
        except:
            continue

    return {'MY_LAYER': dest_id}

The return line is where the crash happens but I'm running out of ideas: by comparison with the template where self.parameterAsSink() returned a tuple (sink, dest_id), I'm guessing the self.parameterAsOutputLayer() returns the dest_id.
Then I try to return the created layer vl and associate it with the dest_id, again, by analogy with the template.

I am stuck here as Qgis (3.0.2) crashes (right at the return line, every line of code gets executed without problem) every time I try to run this current code (but I can make it work fine as a standalone script from the Python Console...)

I have tried the following variations and all lead to invariable crashing:
return 'MY_LAYER'
return {vl: dest_id}
return {dest_id}
return dest_id
return vl
return


texts: QGIS 3 - Consistent symbol size in layout?

texts:
When I'm using the layout tool in QGIS, I would like to be able to preview exactly what the symbols would look like before exporting my map. Currently, symbol sizes are dependent on my zoom scale (i.e. 75%, 120% and so on). The only way that I seem to be able to do this is by setting my symbols to "map units". While a valid option, I don't understand why there isn't a way to preview it within the layout tool? Would it take too long to render?

EDIT: I've included pictures below. The red dots change size depending on scale, but the white squares stay the same because they use map units. The first image is at a 1:1 zoom level:


texts: Inserting string without ' ' in QGIS modeler?

texts:
QGIS 3.4.7

I have a lot of Models where there are some string inputs. In all of them, I have to enter the string between quotation marks ''. I understand the meaning and pourpose of them in any coding.

What I'm looking for is a way that the user doesn´t have to type the '' to avoid mistakes.

How could I concatenate '' with the string typed in the modeler?

That´s how a string input looks in the modeler.

The inputs given in the String Input are used in a simple Field Calculator to fill a new column.


texts: Why does updating Definition Query on layer in ArcGIS 10.1 SP1 no longer allow access to all unique values in dataset?

texts:
In ArcMap 10, setting a definition query on a layer would display only the subset of features:

To add more features to the query would simply involve clicking Get Unique Values to list all the values for that field again:

The steps are the same in 10.1, except when I go to update the query by getting the unique values again, this happens:

The values window no longer shows all the values for that field in the feature class, it only shows the values which were set by the current definition query. I'm not sure if this was a deliberate change (I went from 10.0 SP5 to 10.1 SP1), but it seems that I now have to clear the current query, apply, then go back into the query builder for all values to become available again. Has anyone else encountered this issue?


texts:
I'm working with some Ordnance Survey data, specifically the VectorMap Local (in GML format, October 2018) set and attempting to style the layers using the very handy OS provided stylesheets (here: OS Github Repo link). All the layers render, except for those for which the stylesheets are using rule-based symbology, so this seems to be the source of the trouble. This means that Text.qml, VectorMapPoint.qml, Line.qml all work fine (see repo link above) and that Area.qml, RailCLine.qml, RoadCLine.qml fail to render anything, just white blank canvas.

I've ruled out the obvious candidates. The rule works, or at least testing the rule indicates that it is finding the right features:

I've confirmed that both layer and the whole project are both using the same and right CRS (27700). There aren't any scale related issues (as the screenshot indicates, max/min aren't on very many features, including the one I'm testing). The same layer renders styling Categorised styling using the same column without any issues.

I’ve also tried working off previous commits in GitHub (master, 24 Oct 2018, 12 Feb 2018, 26 Sep 2016) and with other style sets within each of those branches (e.g. “Geopackage stylesheets”), which fails less surprisingly.

Sample data (Vector GML format in my case, os-vectormap-local-gml-sample-data.zip) can be downloaded here: enter link description here


texts: Get minimum, maximum, average, mean with gdal_rasterize?

texts:
I am not 100% sure but it seems like gdal_rasterize takes value of first found feature by default. Is there a way to tell I would like to get minimum / maximum / average / mean (/ count / sum just to get really fancy) of vector features values within cell?

Particularly I have a csv/shapefile  point file with X,Y,Z columns and points were often two or more points fall into one pixel. How can I rasterize maximum value of those within (or minimum, average, mean)?

There are sql "-where" and "-sql select_statement" but not sure if/how these could be used for this case.


texts:
I would like to summarize raster statistics along the edge of each parcel/polygon in a shapefile.

For each parcel shown below in yellow, I would like to calculate the average, min and max value of the raster along the edge bordering the polygon. There will likely be 100,000 polygons that need to be summarized.

Note that the polygons will never overlay the raster, they are only adjacent and share an edge. Ideally those stats could be joined directly to the parcels.

I have ArcGIS 10.2.2 with Spatial Analyst and am comfortable scripting in Python if that is a good path.


texts: Importing XYZ ASCII Grid Data Into QGIS: Why does GDAL resize the raster?

texts:
Question: Why does GDAL import of XYZ gridded data in QGIS modify the dimensions of a raster and can this be avoided?

I'm importing XYZ ASCII gridded data that conforms to GDAL XYZ Gridded Data Standard into QGIS v 2.18.7. The data set represented by the XYZ file has 30 X 30 dimensions where X,Y dimensions of each cell is 0.125 X 0.5 for a total dimension in cells of 240 X 60. X and Y coordinates start at 0.0625 and 0.25, so the outside edges of the cells should span from exactly 0 to 30.

Importing this data through Layer->Add Layer->Add RasterLayer, by drag and drop from the Browser Panel, directly from the OS file browser, or adding via the Georeferencer tool all perform the same task and seemingly import the data properly. However, instead of a 240 X 60 cell raster with 30 X 30 dimensions, these approaches produce a 175 X 175 cell raster with sides measuring 30.012325, and a square cell size of 0.171499.

If I take one of these slightly-off rasters and save it as a GTIFF, the properties within QGIS indicate it is the same, but holding a cursor over the file in an OS file browser window indicates it's 240 X 60, as it should be (see below).

Can someone explain what is going on here and how to avoid it?

Data Example Here

GDALInfo of original ASCII XYZ Grid File: Import_Test_Data.txt

QGIS Properties Metadata of ASCII XYZ Grid File: Import_Test_Data.txt


texts: GDAL does ECW but QGIS (Linux !!) does not load ECW file

texts:
Followed instructions from various posts (each one with slight variations on the theme) and most of the time have no trouble getting the installed GDAL to accept ECW files. However, QGIS stubbornly refuses to accept the fact that GDAL is ECW-enabled and keeps telling me the format is not supported. As no one even mentions the mere fact that QGIS will NOT work even if GDAL itself does, I must be missing something really obvious somewhere (like the need to compile QGIS from source? How to otherwise get ECW in the list of supported raster files in the qgis settings menu...). Grateful if anyone could point that out to me in simple steps, as I've been spending more then a week now on this.

System info:
- Virtual machines using VMWare workstation 7 on Win7 Pro.
- LinuxMint12 (Ubuntu 11 Oneiric) with QGIS 1.7.4 upgraded to QGIS 1.8.0 and GDAL 1.9
- Clean LinuxMint 13 (Ubuntu 12 Precise) with QGIS 1.8 and GDAL 1.9 as installed from qgis package

In all cases gdalinfo --formats list ECW as supported and commandline gdal_translate and gdalwarp work fine with ECW as input and output (translate only), but QGIS knows better and refuses to even try.

Tia, Lewie


texts: Why can not save in opengeosuite ? message: No space left on device

texts:
I have installed opengeosuite in Ubuntu server, and configure to use gwc by leaflet. It was ok and can upload layer using PostGIS, however now it shows this message when try to modify sld style:

  org.apache.wicket.WicketRuntimeException: java.io.IOException: No space left on device

Cannot save anything. Any ideas?

when run in ubuntu: df -h space disk is ok an availble.

Filesystem                    Size  Used Avail Use% Mounted on
udev                          7,9G  4,0K  7,9G   1% /dev
tmpfs                         1,6G  868K  1,6G   1% /run
/dev/mapper/mail10--vg-root   259G   76G  171G  31% /
none                          4,0K     0  4,0K   0% /sys/fs/cgroup
none                          5,0M     0  5,0M   0% /run/lock
none                          7,9G     0  7,9G   0% /run/shm
none                          100M     0  100M   0% /run/user
/dev/sda1                     236M   39M  185M  18% /boot
/dev/mapper/vg_mail-LogVol00  825G  776M  782G   1% /opt


texts:
I am working with multiple grouped layer files that when in the layer file they do not have any of the normal properties. For example, when I click on the properties for an individual layer within the group I only see the "General Tab" but when I drag the layer out of the group I am able to see all of the normal tabs (general, source, selection, etc...).

Anyone have an idea why this is happening?


texts:
I have searched everywhere for an answer to this but as yet have not found one so hoping somebody has done this!
I have a layer of polygons denoting the Parish areas and I have a url in a field that when clicked will open up a web page for that parish and show the carbon footprint e.g. https://impact-tool.org.uk/footprint/footprint?parishId=E04009670&footprintType=territorial&scale=absolute&showSubCategories=true
This is the url in the field.
I want this page to appear in a map tip when hovering over the relevant parish. In the map tip I've used the html notation: <iframe src="[% "comcarboncalcurl" %]" width="640" height="360" border="0"> but it is blank the map tip window shows but is blank.
I've also tried using the url string without reference to the field but again, nothing.
Anyone know the correct syntax or even if it is possible?


texts: Georeferencing xyz layer from Griffiths Valuation of 1850s?

texts:
I have made progress on my project to map land and property owned by individuals in 1850 in Ireland using Griffiths Valuation maps.  I used http://griffiths.askaboutireland.ie/gv4/single_layer/tiles/{z}/{x}/{y}.jpg  with the help of @Steven Kay.  The townlands with beautiful Gaelic names and their borders are shown above in green from OSMIE in comparison with the Griffiths equivalents in red.  There is a difference at the borders which is slightly variable, as expected, from maps separated in time and projection by 160 years.  However, there is a noticeable and fairly constant shift of about 100m. This is most likely caused by something related to projections.

I now intend to create a vector layer which provides the actual property outline of each holder in an overlay.  If I proceed without correcting, the more modern satellite maps like Google Earth and Bing, will not line up properly with the outlines of land and building that I create.

I realize I could use georeferencing to correct this slight misalignment by recording sections of the Griffiths layer as .tif files.  I have already done so with a simple linear transformation of a small area which provides excellent corrected results.  I have tried to insert the corrected projection into the Griffiths layer, but it has no effect.  So, I would like to find a foolproof method of either georeferencing such xyz layers or correcting the projection in the Griffiths layer as it is used.   I just love the convenience of the tiled Griffiths layers I have and the details they provide all over Ireland.

In summary, is there a QGIS process methodology to georeference tiled xyz layers on a server if a final correcting projection is known?

I imagine there is a simple way to do this that a beginner like me is missing.


texts: Categorizing Soil Data using ArcGIS Desktop?

texts:
I have a large shapefile filled with soil information. I am tyring to figure out the best method to classify this soil based on the description into 5 different categories. 1. Clay 2. Sandy 3. Loam 4. Silt 5. Eroded

I have done this manually with smaller datasets. However, I am trying to figure out a way to automate this process.

This is what the dataset looks like:


texts: arcpy export data view to jpeg: label size changes

texts:
I designed a little arcpy script to export different layers in a data view to jpeg (world files included).

It's the first time I actually use arcpy but it works.
The only drawback is that the labels which I set in my mxd, don't come out in the jpeg's as they appeared in the dataview. The labels in my exported jpegs (thus via arcpy) come out very small.

Does anyone know how this can be resolved? I want my jpegs exported with labels, the same way they appeared in the data view. If I adjust the font size in the mxd, nothing happens.

I use Arcmap 10.3 on a Windows 10 laptop.

Below a part of the script which refers to exporting.

if layer.name == "02_GWP_Res":
    layer.visible = True
    arcpy.RefreshActiveView()
    outFile = r"C:\GISlokaal\ArcpyScripts\\" + layer.name + ".jpg"
    arcpy.mapping.ExportToJPEG(mxd, outFile, df, df_export_width=3200, df_export_height=2300, resolution = 150, world_file=True)


texts: Problem when exporting map from printlayout - openlayerplugin disappears

texts: I'm having problems with exporting map from printlayout to PDF or JPG. The Bing map I'm using is disappearing. Can anyone help?


texts: Is there any way to identify user within execute() method of a GeoServer process?

texts: After processing of data within WPS I want to place the result in a folder with a name that corresponds to username, then to send the user URL to access it. To do that, I need a mechanism to determine username within execute() method of my process class. How to do that?


texts: QGIS 2.16 IDs attribute table doesn't match with the row number anymore?

texts: Since QGIS 2.16 the attribute tables count the rows starting with 1, on the previous versions it started with 0, so the $ID could match with it. Is this some kind of mistake? Do I have to use another command? Don't really understand the usefulness of this.


texts: Leaflet.TextPath: how to set zoom level from which text would be displayed

texts:
I am using the Leaflet.TextPath plugin for labeling polylines. It works fine into my Leaflet map, but I want to set the zoom level (scale-based visibility) from which labels would be displayed on the map.

How can I do that?


texts: Getting unicode encode error in ArcPy/Python script?

texts:
When I am running a script (it checks to which database a mxd is connected) I get the following error:

return codecs.charmap_encode(input,errors,encoding_table)
UnicodeEncodeError: 'charmap' codec can't encode character u'\u2265' in position##: character maps to <undefined>

Script:

import os, arcpy

directory = r"C:\Project"
for root, dirs, files in os.walk(directory):
for myFile in files:
        fileExt = os.path.splitext(myFile)[1]
        if (fileExt == ".mxd"):
            fullPath = os.path.join(root, myFile)
            print myFile

            myMap = arcpy.mapping.MapDocument(fullPath)
            layers = arcpy.mapping.ListLayers(myMap)

            for layer in layers:
                wsPath = ""
                try:
                    wsPath = layer.workspacePath
                except:
                    pass
                if wsPath:
                    print "    " + layer.name + " : " + wsPath

Is there a way to keep the script running and catch an exception?


texts: Custom Tool for creating a cirle out of straight segments

texts: I used to use an old tool coded by somebody that would create a circle much like the current ESRI Circle tool, but it would make the circle out of straight segments with many vertices.  The dll can not longer be registered as I believe Arc has just outgrown it (tool was made in 2003). Is anybody aware of another tool that can accomplish this same feat? The arc segments generated by the ESRI circle tool cause headaches with some of our Topology checks.


texts: GIS or Remote Sensing technique to map sub-surface features?

texts:
I am trying to delineate agriculture tile drainage in Midwest.

Does anyone know of any GIS and/or sub-surface remote sensing techniques to detect tiles/objects buried in shallow depth of soil?

Does anyone have any suggestions as to what type of remote sensing can be used to detect subsurface features?

I can think of ground penetrating radar (GPR).


texts: Exporting raster data to a spatial grid in R

texts:
I have a raster as follows:

library(geostatsp)
library(raster)
library(sp)

xmin <- -11.69
xmax <- 2988.31
ymin <- -749.97
ymax <- 1650.03
res_r <- 50
r <- squareRaster(extent(xmin, xmax, ymin, ymax), cells=xmax*ymax)
res(r) <- res_r
values(r) <- runif(ncell(r))
print(r)

class      : RasterLayer
dimensions : 48, 60, 2880  (nrow, ncol, ncell)
resolution : 50, 50  (x, y)
extent     : -11.69, 2988.31, -749.9699, 1650.03  (xmin, xmax, ymin, ymax)
crs        : NA
source     : memory
names      : layer
values     : 0.0005739587, 0.9996859  (min, max)

Using the raster, I can obtain the cell center coordinates as follows:

coord_r <- rasterToPoints(r, spatial=TRUE)
map_pts <- data.frame(x=coordinates(coord_r)[,1], y=coordinates(coord_r)[,2], coord_r@data)

However, I would like to obtain the coordinates of cell vertices (or four corner coordinates of each grid cell) instead of cell center coordinates. To do this, I've created a spatial grid by offset, cell size and dimensions as follows:

r_g <- GridTopology(c(xmin,ymin), c(res_r, res_r), c(ncol(r) + 1, nrow(r) + 1))

The grid r_g should contain 2989 observations.

But I don't know how to export raster data to the grid r_g. Any help would be greatly appreciated

Edit:
Here are some details about how I obtain a grid with 2989 observations:

## 1) I retrieve cell vertices
polyg <- rasterToPolygons(r)
polyg@data$ID_poly <- rownames(polyg@data)
colnames(polyg@data) <- c("variable", "ID_poly")
polyg_vertices <- spbabel::sptable(polyg)
polyg_vertices <- polyg_vertices[!(polyg_vertices$order_==5), c("object_", "order_", "x_", "y_")]
colnames(polyg_vertices) <- c("ID_poly", "order", "x", "y")
tab_polyg_vertices <- join(polyg_vertices, polyg@data, by="ID_poly")

## 2) I delete duplicated values
tab_polyg_vertices <- tab_polyg_vertices[,c("x","y")]
tab_polyg_vertices <- unique(tab_polyg_vertices)
dim(tab_polyg_vertices)
[1] 2989    2

However, I don't know how to retrieve raster data (here, the values of the column "variable") because there are several values for the same corner coordinate.


texts: ArcGIS - How to programmatically create new dynamic layer for service?

texts:
I'm trying to create a number of widgets that will visualise different data by geographic area, e.g., the Potential Losses for Instrument X in Year Y for the Counties of State Z.

Ideally this would be a map of the state (tick can do that)
with the showing the county boarders (use javascript api to display basemap layer )
and then colour each county according to loss. :(

Going down the ArcGIS route I am assuming I need to create a dynamic layer here with all of the polygons for each county of each state defined and then somehow colour them according to data pulled out of a DB showing loss for that instrument, for that county for that year.

Could someone please

A) validate that that is the correct way to go and
B) point me in the right direction to some resources to help me learn
how to achieve this.


texts: How to force GeoServer to remove, or change, NoData Value from GeoTiff?

texts:
I am trying to display vector data of state borders, downloaded from GeoServer as GeoTiff, in QGIS.
By default downloaded raster looks like this:

After disabling NoData value in raster->Properties->(uncheck)NoData Value:

Current SDL code is:

        <FeatureTypeStyle>
            <Rule>
                <Name>Ime</Name>
                <MinScaleDenominator>70000</MinScaleDenominator>
                <MaxScaleDenominator>170000</MaxScaleDenominator>
                <PolygonSymbolizer>
                    <Fill>
                        <CssParameter name="fill">#FFFFFF</CssParameter>
                        <CssParameter name="fill-opacity">0.1</CssParameter>
                    </Fill>
                    <Stroke>
                        <CssParameter name="stroke">#000000</CssParameter>
                        <CssParameter name="stroke-width">0</CssParameter>
                    </Stroke>
                </PolygonSymbolizer>
                <TextSymbolizer>
                    <Geometry>
                        <ogc:Function name="centroid">
                            <ogc:PropertyName>SHAPE</ogc:PropertyName>
                        </ogc:Function>
                    </Geometry>
                    <Label>
                        <ogc:PropertyName>NAZIV</ogc:PropertyName>
                    </Label>
                    <Font>
                        <CssParameter name="font-family">Verdana, Geneva, sans-serif</CssParameter>
                        <CssParameter name="font-size">16</CssParameter>
                        <CssParameter name="font-weight">bold</CssParameter>
                    </Font>
                    <LabelPlacement>
                        <PointSymbolizer>
                            <AnchorPoint>
                                <AnchorPointX>0</AnchorPointX>
                                <AnchorPointY>0</AnchorPointY>
                            </AnchorPoint>
                        </PointSymbolizer>
                    </LabelPlacement>
                    <VendorOption name="autoWrap">80</VendorOption>
                    <VendorOption name="maxDisplacement">20</VendorOption>
                </TextSymbolizer>
            </Rule>
        </FeatureTypeStyle>

Is there any way to disable NoData Value in Geoserver?


texts:
I am looking for regionalization software, that is, I want to cluster regions based on multivariate attribute data with the added constraint of spatial contiguity.

I have found various software for doing this including those listed on another GIS.se question

All so far seem to standardize data before clustering, so data measured in different units are treated equally. None so far, however, offer the ability to weight variables. So for example if I determine that population density is twice as important as house price in determining regions, I would like to be able to specify this to the algorithm.

Is there some regionalization software out there that allows variable weighting?


texts: What's the best way to store State Plane data in SQL 2008R2?

texts:
I have a client with data that is currently in a State Plane coordinate system.  They are looking to store this data in SQL Server 2008R2.

SQL does not "support" the SRID for the State Plane in question.

I've seen some people suggesting using geometry instead of geography for storing State Plane data.

I also know that I could convert the State Plane data over to WGS84 if necessary.

What I'm wondering is what's the "best" way to store State Plane data in SQL 2008R2?


texts: I've got a polygon shape layer, and each polygon has 200 or more vertices. I need to divide every polygon by attributes (vertices). Each polygon can not have more than 80 vertexes. I only found plugins that divide polygon by area.


texts: Trend (Spatial Analyst) - Chi-square value

texts:
I have to calculate the slope of a line based on a dataset of points with a height value.

I am doing this by calculating the slope over a raster created using the trend function (linear 1st order). when i create the raster using the trend function I get also an RMSE value and the chi^2 value. I understand how the RMSE value is created, but I dont understand how the chi^2 is calculated.
Is there someone who can explain to me how this value is calculated, and/or can tell me how to find out the degrees of freedom so i can interpret this value?


texts: Mapserver7 PostGIS itemnquery problem migrating from 6.2.1

texts:
I am migrating an existing application from Mapserver 6.2.1 to Mapserver 7.0. How did the expression-Querying change there?

When doing an itemnquery or itemquery it always fails in querying my postgis-Layers.
This is my query:

mapserv.exe "QUERY_STRING=mode=itemnquery&schluesselgesamt='0335502400092'&qstring=hausnummer='5'&qlayer=alkisadresssuche&mapsize=1129+646&map=E:\tmp\tmp_alkis.map&"

I didn't use any FILTER in my layer but it seems that it therefore includes an empty filter which can be found in the Postgres-log:

syntax error at or near ")"

here is the sent query from Postgres-Log

select "schluesselgesamt","id","hausnummer",encode(ST_AsBinary(ST_Force2D("geom"),'NDR'),'hex') as geom,"id"
from (SELECT schluesselgesamt, id, geom geom FROM alkis.adresssuche ) as ergebnis where geom
&& ST_GeomFromText('POLYGON((593025 5891783,593025 5893724,595424 5893724,595424 5891783,593025 5891783))',25832)
and () and (schluesselgesamt='0335502400092')

It works with ONE query parameter but not with multiple parameters. I tried several different calls to mapserver and it works without EXPRESSION or FILTER or PROCESSING=NATIVE_FILTER.

This was the important part in the original query: :&qstring=schluesselgesamt='0335502400092'+AND+hausnummer='5'&

The Mapserver homepage just told me, that the FILTER has changed but not that it is now needed. Also is there no hint on the Mapserver EXPRESSION page on possible changes. The Mapserver7 changelog only gives some tickets about changes in query but they didn't help much (http://mapserver.org/development/changelog/changelog-7-0.html#changelog-7-0-0-beta2)

Does somebody has any hints for me in this case? Any not so well documented changes in Mapserver7?

By the way: here is my layer:

LAYER
    CONNECTIONTYPE postgis
    NAME "alkisadresssuche"
    CONNECTION "user=***** dbname=local_database host=localhost password=****"
    DATA "geom from (SELECT schluesselgesamt, id, hausnummer, geom FROM alkis.adresssuche
            ) as ergebnis using unique id USING srid=25832"
    STATUS On
    TYPE POLYGON
    DEBUG On
    HEADER   "../templates/utf8/search_h.html"
    TEMPLATE "../templates/utf8/search_q.html"
    CLASS
       NAME "Alkis-Adresssuche"
       STYLE
          WIDTH 5
          OUTLINECOLOR 0 0 255
       END
    END
    VALIDATION
        "schluesselgesamt" "."
        "hausnummer" "."
        "qstring"  "."
    END
END


texts: Eps files from map export are not containing a vector image

texts:
I have made a custom map export service, which exports my map in a range of formats.
However I have noticed that .eps files, exported using the arcPy function:

  arcpy.mapping.ExportToEPS(mxd, Output_File, resolution=96, picture_symbol='VECTORIZE_BITMAP', convert_markers=True)

are not a vector images. The map exported, is not well scalable and it seems pixelated.  Also, we can see that tiled layers and basemaps are not exported as vector images. However in feature/ graphic layers that use picture markers and picture fills will be converted to vector. Any ideas why is that happening? Tile layers cannot be exported as vector images?


texts:
I am trying to use the GetFeatureInfo functionality of Openlayers v3.

I am taking the layer name dynamically from the layerswitcher and trying to  pass it into WMS like:

var temp=JSON.stringify(LayerName);
var newTemp = temp.replace(/"/g, "'");
var wmsSource = new ol.source.TileWMS({
      url: 'http://localhost:8080/geoserver/wms',
      params: {'LAYERS': newTemp},
      serverType: 'geoserver'
    });

Is it possible to pass a JavaScript object to "Layers:"?

If I pass the layer name directly it works smoothly but I need a dynamic functionality.


texts: PyQGIS 3: How can pie charts have a fixed size when automating the process?

texts:
I'm unsure how to automate the process for PyQGIS to update the pie charts to a fixed size of 12.
The approach I used is not working out and no pie charts are being created.

lyr = iface.activeLayer()

diagram = QgsPieDiagram()

ds = QgsDiagramSettings()

dColors = {'Field1': QColor('#5a09a6'),'Field2': QColor('#ff9200')}

ds.categoryColors = dColors.values()
ds.categoryAttributes = dColors.keys()
ds.sizeType = 0
ds.categoryLabels = ds.categoryAttributes

dr = QgsLinearlyInterpolatedDiagramRenderer()
#dr.setClassificationField('SAMPLE NUM')
dr.setDiagram(diagram)
dr.setDiagramSettings(ds)

lyr.setDiagramRenderer(dr)
dls = QgsDiagramLayerSettings()

lyr.setDiagramLayerSettings(dls)
lyr.triggerRepaint()


texts: Invalid QgsRasterLayer when attempting to load a Google Static Maps API image into QGIS

texts:
In the QGIS Python console, this block of code returned a valid object in QGIS 2.16.3 and I think 2.18.6. I could view the static image in QGIS.

>>> reg = QgsMapLayerRegistry.instance()
>>> raster = QgsRasterLayer('https://maps.googleapis.com/maps/api/staticmap?center=51,0&zoom=10&size=2550x3300', 'test')
>>> raster.isValid()
>>> reg.addMapLayer(raster)

At least since 2.18.9 (tested 2.18.15 as well), raster.isValid() returns False.

This is not specific to the Google Static Maps API. The code block also used to work if loading other PNG/JPG images.

Am I missing a required component?


texts: Label sums of surfaces grouped by common attribute values

texts:
I want to label the sums of surfaces ("Groesse") grouped by common attribute values ("Block") in QGIS 3.16. I need one label for each group.
The following pictures help to describe the expected result.
Here is an example of my attribute table:

Here is the result I want to get. 5 sums grouped by the "Block" column are represented below: 3 orange (D1, D2 and D3), 1 grey (S1) and 1 green (G1)

I thought to create a virtual layer with a function that sum the surfaces ("Groesse") per attribute value("Block"). Despite my attempts, I do not know how to write it correctly in SQL (sum("Groesse","Block") in Python).

I also tried to display labels through a Rule-based expression in using the aggregate function, without success.


texts: Empty result layer was created by raster calculator

texts:
I'm a new user of QGIS and learning some tips from this book https://www.amazon.com.au/Learn-QGIS-Andrew-Cutts/dp/1788997425. I'm currently practising how to use raster calculator to find areas at elevations > 1000 m with a slope < 5 m.

First, I downloaded the image from http://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/TIFF/srtm_05_01.zip (WGS84 projection)
Then I loaded this .tiff file to QGIS and saved a new layer DEM_local with another projection (NAD27-ESPG:2964)

After that I used Slope tool for this DEM_local to make a new slope layer.

Use raster calculator to find elevations > 1000m and Slope < 5m

Although I followed book guides the result layer was still disable.
Some previous posts said it's caused by improper layer names or projections. As they suggested, I've changed srtm_05_01 to srtm0501 or checked projections , but the problem still remains.

The range of raster_calculator_eg should be 0 - 1 (similar to this book) instead of +/- 1.79769e+308

Could anyone please help me?


texts:
All my layer URLs and Feature Layers are stored in an array.

with an foreach loop, I am creating the FeatureLayers from the URLs,

dojo.forEach(window.layerUrls, function (info, idx) {
            window.featureLayers[idx] = new esri.layers.FeatureLayer(
            window.layerUrls[idx], {
                infoTemplate: infoTemplate,
                mode: esri.layers.FeatureLayer.MODE_ONDEMAND,
                outFields: window.outFields
            }
          );
            map.addLayers([window.featureLayers[idx]]);
        });

This is the code to generate the legend:

dojo.connect(window.map, 'onLayersAddResult', function (results) {
            var legendDiv = dijit.byId('legendDiv');
            if (legendDiv) {
                legendDiv.destroyRecursive(true);
            }
            var layerInfo = dojo.map(results, function (layer, index) {
                return { layer: layer.layer, title: layer.layer.name };
            });

            var legendDijit = new esri.dijit.Legend({
                map: map,
                layerInfos: layerInfo
            }, 'legendDiv');
            legendDijit.startup();

        });

Whe I run the page this is what the legend looks like and I have no idea why:
http://abload.de/img/legend22uj3.jpg
Here is the code related to this picture:

Keine Legende vorhanden

Legende wird erstellt...

Legende wird erstellt...

Legende wird erstellt...

Legende wird erstellt...

Legende wird erstellt...

Additonally: when zooming the legend disappears and only one random FeatureLayer-legend remains

Where did I go wrong and why?