texts: Filter (query) OSM .pbf file in QGIS graphical modeler or Python (QuickOSM or related tools)

texts:
I would like to automate a simple operation in the QGIS graphical modeler or a simple Python script.

add all line features from a .pbf file
filter all features with certain values within a certain field (e.g. "waterway" = 'canal')

Doing this manually is fairly easy, by just dragging the file in the canvas, selecting all line features, and then right clicking it to select filter.
However, in the graphical modeller I struggle with the very first step of loading the .pbf files' lines.
QuickOSM
I know there is QuickOSM, but again I can do it manually, but the algorithm "Open sublayers from an OSM file" that you can add to the graphical modeler doesn't do the job. First you cannot distinguish between lines and other features types (as in the manual tool). And then, I can link it to some filter algorithms in the graphical modeler (eg. "Filter by geometry type", "Extract by attribute", or "extract by expresion"), but the best I get is an empty line feature :-(.
It is indicating separate outputs for the different feature classes, but how can I access them?:
It seems that it is possible filtering the pbf file while loading it by configuring the .ini file (https://gdal.org/drivers/vector/osm.html). Thereby, I hope to obtain only one line feature output that I can use as input for follow up algorithms, but no idea how to do it.
esy-osmfilter
Besides, esy-osmfilter is an easy Python solution, but it runs into memory crash for large data extracts.


texts: Is there an issue with Zonal Statistics calculating wrong mean values?

texts:
I tried to solve the following issue by myself but without success. On the one hand I have a raster file which is the result of the processing "shapes to grid" tool. On the other hand I have a boundary vector file. Now I wanted to use the zonal statistic tool via processing toolbox. But the result isnt correct. The mean value is not what it should be. The metadata (grid file proberties) shows me the correct mean value. The crazy thing is, when I convert the raster into a shape file again and then converting this shape into a raster file once more the zonal statistics tool works fine.

Could please anyone check this for me? I'm using QGIS 2.2. The files are here (EPSG 25833):
Link


texts:
What is the function for getting SRID from Geometry field.

I have inserted a shapefile into SQL Server table with a Geometry field defined. I know the projection of shapefile is Swreff99_1800 (EPSG:3011).

Now my question is How I get this SRID:3011 by using some OGC method of the geometry field?


texts: Creating geodatabase, that has all layers used in creating map from different geodatabses?

texts:
I have been working on creating a city map, that had layers sourced to different gdbs stored on my desktop, now my final product is ready after 3 months and I have to package the datasets and send only one Geodatabase that consists of all the layers used to create the exact city map. How can I do that without disrupting the original layout of the map?

Pictures for reference. (Sorry, but my boss wouldn't let me reveal the gdb names)


texts:
This seems too simple to be asking this but I am really new to PostGIS/Postgres have been going around in circles with this. I am performing a really simple spatial query where I would like to make a subset of a table of points conditional on whether they are contained within a polygon in another table.

I can return the rows I want from my table, however, they are being joined to the other table. All I want is a subset of one table based on the condition, not a join per se.

Here's my query:

SELECT *
FROM
public."KDhh_survey",
public."ur_pilot_survey"
WHERE
ST_contains(public."ur_pilot_survey".the_geom, public."KDhh_survey".the_geom);

I only want the records in table public."KDhh_survey" returned in a new table or view, but ALL records are returned from both and a join created.

I promise I have searched high and low for an answer to this, but every example I found for ST_contains is based on spatially joining two table together. I really only want to query my data set, so that I can create a subset of my original table.


texts: How to add randomly distributed points within boundaries of zip codes?

texts: I’m creating a map to show the distribution of students in my school’s department. I have a .shp file from the US Census displaying zip code areas in the vector layer. Now I have an Excel spreadsheet that lists how many students are in each zip code per academic major. I was wondering if there was a straight forward way to create a vector layer displaying the appropriate number of points evenly distributed throughout each zip code. I tried to join the student data converted to a .cvs file to the Census zip code vector layer. But nothing appeared to have happened. Thanks.


texts: How to remove duplicate points using ArcGIS ModelBuilder?

texts:
I want to remove duplicate points and I have this model.

I want to deselect all feature that have the same value of the OID. So I want to add value of the OID in the query of Select by Attribute. How do I incorporate this value in a SQL query ???

I try : "%Feature OID%" = "OID"   (does not work)
and a lot of other tests and nothing work

Someone have an idea?


texts: Calculating distance in m in xyz between GPS coordinates that are close together

texts:
I have a set of GPS Coordinates and I want to find the speed required for a UAV to travel between them. Doing this by calculating distance in x y z and then dividing by time to travel - m/s.

I know the great circle distance but I assume this will be incorrect since they are all relatively close together?

Is there an accurate way to do this?


texts: Join attributes from multiple polygons to an overlaying polygon

texts:
I have a polygon-shapefile, containing 562 single features (feeding -layer). I want to add their attributes to an overlapping polygon-layer (bear_buffer -layer), containing 1330 features. The bear-buffer features are overlapping in most cases, so I decided to use the intersect-tool for both layers, because I also need the exact area-sizes of the "feeding"-features inside the single "bear_buffer"-features (always total). This resulted in about 7,030,000 single features.

feeding -layer

bear_buffer -layer

both layers

Now I want to add several attributes of the "feeding"-layer to the "bear_buffer"-layer, so I get for every "bear_buffer"-feature the attributes of the "feeding"-layer. The following 2 pictures will give an example of the 2-input-tables:

example attributes feeding -layer

example attributes bear_buffer -layer

In the end, the result should be like on the following picture, so that I can import the table easily into Excel (because this doesn't work with the 7 millions features of the intersect-ouput) and doing my analyses.

example output-table (AAA, BBB, ... just represent other attributes), "bear_buffer"-features will be repeated for every "feeding"-features joined to the table

So the question is: How can I do this? Have tried a lot of opportunities, nothing helped me to achieve the perfect result. If Excel could handle with 7 million features (rows) it would be easy with pivot table, but in this cases it's a bit more difficult.


texts:
Is it possible to connect to a SQL Server database with GRASS?

I tried this way:

>db.connect driver=ODBC database='dbname=mydb,host=myserv'
(Mon Dec 12 15:34:27 2016) Completed (0 sec)`

But next, when I try to show a list of tables:

>db.tables
-p
ERROR: Unable to open database <dbname=mydb,host=myserv>

The database is hosted on a server on my local network and I use GRASS 7.0.5 and SQL Server 12


texts: how we create line between two cordinates on polymap?

texts: I am working on d3 and polymaps, I plot map and create circle on the given latitude and longitude which retrieves from the JSON file but now I want to connect those circles with line. I tried to do it with SVG but I am unable to find the accurate result.


texts:
Is there a programmatic way to clear the cache of a map in ArcGIS Engine 10.0?

It is apparent that there is some caching scheme operating under the covers, and I'd like to clear all of its cached data.  To be clear it isn't leaking memory, I just don't want it to be caching anything.

This is a fairly vanilla (*.mxd).  Just a raster dataset, no basemap layers or anything fancy going on.


texts: how to Replace masked out cloud areas with another satellite image?

texts:
I have a Sentinel 2 image which has heavy cloud cover. I would like to mask out the clouds, and replace these areas with a Landsat 8 image.

(1) How do I mask out cloud?

(2) how to I replace the cloud areas with the Landsat 8 image?

I am using SNAP software.


texts: Rasterizing Single Shapefile Based On Its Attributes

texts:
I am trying to rasterize single shapefile based on its attributes automaticly. Let us suppose that I have shapefile at least 50 attributes. I would like to rasterize this shapefile based on its unigue field attribute in automatic way. It means that every unigue attribute will be stored in single raster band with values of pixels of this attribute. I tried to write an R script, but I was not successfull. Here is part of my code:

library(raster)
library(rgdal)

# setting working directory
vstupni_adresar <- choose.dir(getwd(), "Zadej cestu k pracovnímu adresáři")
prac_adresar <- setwd(vstupni_adresar)

# loading raster data
r <- brick(choose.files(getwd(), "Vyber obrazová data"))

# loading vector data
shp <- shapefile(choose.files(getwd(), "Vyber polygonová data"))

# resolution and extent of raster
cellsize <- 20
ext <- extent(r)
rs <- raster(ext, res=cellsize)

# loop of rasterization
for (i in 3:names(shp)){

# vlastni proces rasterizace
rasterizace <- rasterize(shp, field=i, rs)
plot(rasterizace)

# name of output raster
outname <- paste("B", i, ".tif", sep="_")

# export raster
rase <- writeRaster(rasterizace, filename=outname, format='GTiff',
                    overwrite=T)
}

But script produces an error. For the matter of speeding up computations I would like to add parallel implementation with foreach package as well.

As for the task itself I does not have to be R, I tried to do this with gdal_rasterize utility which is really fast, but I did not succeed. I prefer R to anything else because it is the thing of which I am the most familiar with.

In another way. For example when I want to rasterize polygon shapefile with attribute called ,,AREA" and save it as a raster file named ,,RASTER" it will look like this:

gdal_rasterize -a AREA -tr 20 20 -of Gtiff shapefile.shp RASTER.tif

But there is a situation when one would like to rasterize the whole attribute table it means to create separate command or batch file with commands above where only parameter -a will differ for the whole attribute table based on name of attribute column. Simplier and much more effective way to do so is to create a script which will loop through list of all present atrributes of polygon shapefile.


texts: Using GDAL/Python to stack georeferenced images of different sizes

texts:
I am in the process of porting a code I wrote in IDL (interactive data language) to Python but am running into a bit of a problem that I am hoping someone can help me with.

The code goes like this:

take individual classified Landsat GeoTIFFs (say there are N individual 1-band files per scene, each representing a different day) and further reduce these images to three binary-themed 1-band images (water and not water, land and not land, water/land and not water/land). This will be done by reading the rasters as matrices and replacing values.
** I don't actually need to have these images, so I can save them as memory or just keep them as numpy ndarrays to move to the next step
stack these images/arrays to produce 3 different (1 for each 'element') N-band stacks (or a 3-dimensional array-- (samples, lines, N)) for each scene
total the stacks to get a total number of water/land/water&land observations per pixel (produces one 1-band total image for each scene)
other stuff

The problem I am running into is when I get to the stacking, as the individual images for each scene vary in size, although they mostly overlap with each other. I originally used an ENVI layer-stacking routine that takes the N different-sized 1-band images for each scene and stacks them into an N-band image with an extent that encompasses all of the images' extents, and then reading the resulting rasters in as 3-D arrays to do the totals. I would like to do something similar with GDAL/Python but am not sure how to go about doing so. I was thinking I would implement GDAL capabilities of GeoTIFFs by using the geotransform info of the images to somehow find the inclusive extent, possibly padding the edges of the images with 0's so they are all the same size, stacking these images/3-d arrays so that they are correctly aligned, then computing the totals. Hopefully there is an easier way, as I'm not sure how to pull that off.

Does anyone have any suggestions or ideas as to what would be the most efficient way (or, any way really), to do what I need to do? I'm open to anything.


texts: PyQGIS spatial query geometry function disjoint incorrect output

texts:
legend=iface.legendInterface()
layers=legend.layers()
layer1=layers[0]
layer2=layers[2]
pointCounter=0
for pointFeat in layer1.getFeatures():
    pointGeom=pointFeat.geometry()
    for polyFeat in layer2.getFeatures():
       polyGeom=polyFeat.geometry()
       if pointGeom.disjoint(polyGeom):
          layers[0].select(pointGeom.id())
          pointCounter=pointCounter+1

print pointCounter

not getting the correct output for disjoint


texts: Vector layer no longer editable in QGIS?

texts:
This problem was previously mentioned for QGIS 2.16 but I was hoping that QGIS 3.16.5 will not have the same problem. Somehow a vector shapefile gets into a state where it is not possible to edit the layer (the pencil is grey in the attribute table).
I was adding lines to a vector layer and saved the changes but nothing was saved (all changes were lost) and now the layer is not editable. Restarting QGIS does not help. Removing the layer and adding it again does not help.
I previously I had the same problem and I had to create a new vector layer.
Is it a bug or a mistake on my side?
The problem shapefile can be downloaded from https://www.dropbox.com/s/2w52aurd9tbi2hr/Uneditable%20layer.zip?dl=0


texts: Open ArcGIS Geodatabase as python file object

texts:
I'm developing a script to POST a .GDB to a geoprocessing service. Part of my workflow involves opening it while building the POST request:

files = {'file': (open(C:\\Temp.gdb, 'r'))}
r = requests.post(uploadsUrl, files=files, data=data)

I'm able to execute the first line with a personal geodatabase (.mdb) but receive an error when opening a GDB:

IOError: [Errno 13] Permission denied: 'C:\\Temp.gdb'

Is it possible to use the open() command on a file geodatabase? I could read through all feature classes and POST them individually to the REST API, but moving the GDB would be much easier.


texts: Where to find Data with site that have hazardous, toxic or radioactive substances?

texts:
Does anybody know where to find data that lists the location of sites that have hazardous, toxic or radioactive substances?

I'm working on a HUD project and need to be able to find out if sites for new housing is within a certain distance of gas stations or other places that are potentially hazardous.

Specifically interested in Florida, USA.


texts: Filter QgsAttributeTableModel without geometry

texts:
Is it possible to filter a QgsAttributeTableModel without geometries.
I was not able to filter a QgsAttributeTableModel without geometries, so I added a fake geometry column to my table ("geom_fake"). Now it works, as I want.

I am working on a QGIS plugin and I want to display data in tables, sometime with and sometime without geometries. So I use a QgsAttributeTableModel to achieve this.

Is there a smarter solution to filter tables without geometries (like in QGIS, where you can filter tables with or without geometries very easy).

Addition:

I am using QGIS 3.4 with a SpatiaLite-database (SQLite) created with QGIS.
Here is a summary of the used code (this is not a working example).

first i set the model and the tableview:

    self.layer_uri = '''dbname='C:/Daten/db/GisDaten.sqlite' table="alm_gdb_main" (geom_fake) sql='''
    self.layer = QgsVectorLayer(self.layer_uri, 'GDB-Grundstücke', 'spatialite')
    self.layer.getFeatures()

    self.canvas = iface.mapCanvas()

    self.vector_layer_cache = QgsVectorLayerCache(self.layer, 10000)
    self.attribute_table_model = QgsAttributeTableModel(self.vector_layer_cache, self)
    self.attribute_table_model.loadLayer()

    self.filter_mod = QgsAttributeTableFilterModel(self.canvas, self.main_table_model, self)

    self.uiMainListTblvw = QgsAttributeTableView(self)
    self.uiMainListTblvw.setModel(self.filter_mod)

to filter the model call the following code:

    selection = self.layer.getFeatures(gsFeatureRequest().setFilterExpression(self.filterString.text()))
    ids = [i.id() for i in selection]
    self.filter_mod.setFilterMode(3)
    self.filter_mod.setFilteredFeatures(ids)

'self.filterString' is a lineedit where you can enter the filterstring for testing.

This code above works as i want, when the layer has a geometry column.
To handle a layer without a geometry-column i would change 'self.layer_uri' to '''dbname='C:/Daten/db/GisDaten.sqlite' table="alm_gdb_main" sql='''.
But now the filter doesn't work correctly anymore.


texts:
I'm trying to modify the scale of the map so I can further export to pdf.

The mxd document has a fixed map scale option set. The initial value is 250 and final is still 250 unchanged. The code looks like this:

mxdPath = GetResPath() + "\\Template.mxd"
mapDoc = NewObj(esriCarto.MapDocument, esriCarto.IMapDocument)
mapDoc.Open(mxdPath)

pActiveView = mapDoc.ActiveView
pFeatureLayer = pGPUtil.OpenFeatureLayerFromString(path)
pLayer = CType(pFeatureLayer, esriCarto.ILayer)

pMap = mapDoc.Map[0]
pMap.MapScale = 500.0 #this will change the value, but the mxd just doesn't care
pMap.AddLayer(pLayer)
pActiveView.Refresh()
outDocumentPath = GetTempPath() + "\\Document.mxd"
mapDoc.SaveAs(outDocumentPath)
mapDoc.Close()

I also want to mention that setting the Extent through pActiveView.Extent works, but it's harder to get an extent out of a scale

EDIT 1: I have already tried with both Automatic and fixed, on Automatic if I try to get the property it will throw and Exception and if I try to set the property nothing happens also.

EDIT 2: It seems the ArcObjects does actually modify the MapScale when setting the property, but it won't save this option when calling mapDoc.Close(). The workaround was to load the .mxd again using ArcPy and setting the .activeDataFrame.scale property.

I hope this is helpful for people having related problems!


texts: How to split lines where a different set of lines intersect them, no vertices

texts:
How do I use one shapefile to split the lines into another shapefile?

I have a line shapefile of streams and another line shapefile of roads. I would like to split the streamlines where the roads intersect with the streams. These two shapefiles don't share any vertices.


texts: Searching locations in any ArcGIS Basemaps?

texts:
I am wondering if it is possible to search locations in the basemap that is used as a background map in my project or whether I can add any searchable layer (such as it happens in Google Maps?).

The problem is that I need to edit a set of points that are currently located in the wrong places. As I don’t have an extensive knowledge of the region under study, and it would be very cumbersome to first search the location in Google Maps, and then go back to ArcGIS and move the points.

Does ArcGIS have any search tool that finds the addresses/placenames of the locations I am looking for in the basemaps.


texts: Count points in polygon grouping by attribute using PyQGIS

texts:
I have a point layer with an attribute table like this:

id_point  class
       1      A
       2      B
       3      A
       4      C
       5      B
...

What I need is to count how many points of any class are in every polygon. Something like this

id_pol   A   B   C
  pol1   1   2   0
  pol2   1   0   1
...

Is there a way to count points within polygons grouping by a point attribute using Python in QGIS 3?

I think QGIS 2 had a Count Unique Features by Polygon but I can't find it in QGIS 3.


texts:
I have been trying to convert a float value, for instance, 6.764789 minutes, to a time formated number 00:00:00.

First, I created a New Field of Time, and then used this two ways (using the Field Calculator):

def FormatTime(st_t):
    sec = st_t*60
    m,s = divmod(sec,60)
    h,m = divmod(m,60)
    return "%d:%02d:%02d" % (h, m, s)

And just making a simple operation:

New Time Field = Float Number Field/ 60

Both are correct, but results are 00:00:00


texts:
I have an infrastructure layer that is lat/long values and has an SRID of 4269.

I have brought this shapefile into PostGIS to try and reproject the data to 3424 (state plane of NJ)
this is my query to reproject

ALTER TABLE infrastructure
   ALTER COLUMN geom
   TYPE Geometry(Point, 3424)
   USING ST_Transform(geom, 3424);

but it still shows up as lat long values when I do the below query

I then brought the infrastructure shapefile into qgis and all the hundreds out points showed up in one area outside of NJ...I am not sure what i am doing wrong here


texts: Is it possible to indicate an angle to a point marker symbol from field without using representation?

texts: Is it possible to indicate an angle to a point marker symbol automatically from data from a field without using Representation ?


texts: Geoserver SLD stylesheets - what does "se" namespace mean?

texts:
I frequently use QGIS to generate SLD stylesheets for geoserver layers. When I do it that way, namespace with value "se" is added to SLD tags.

Example:

<se:Rule>

When I make a stylesheet manually, I don't use namespace, and it works the same. Why is namespace added in this case?


texts: How to merge overlap polygon respectively in one shapefile in QGIS

texts:
I want to merge the polygons of one shapefile if they overlap， but I don't know what to do.

I try'dissolve' but get only one final polygon. This is not what I need. I want to merge them respectively.


texts: How to move polygon fit to Google satellite in Qgis

texts: I have a huge polygon layer of land areas and houses in a river delta. It does not really align with the images form Google satellite, so I tried moving the polygon layer. When I do so, it can't be loaded.


texts: Why are my symbols not displayed correctly in the legend and after a PDF-Export?

texts:
I work with QGIS 2.0.1.

I have two problems with the representation of symbols.

The first problem is that after I selected a symbol with the type: font and adjust the size of the symbol with the unit: map units, the symbol is not displayed completely within the Legend.

Is this a general problem that the icon size is also scaled within the legend?

Legend:

The next problem is that when I want to export my map to a PDF file the Symbols are not displayed correctly.

The icon consists of two components.

I have tried two different versions but the result was not useful.

Here are the two variants i have tried:

Symbol as font and white background as a Font.
Symbol as font and white background as a simple symbol with a square.

In QGIS both variants provide a correct representation.

But the result was similar for both variants. The Symbol components are shifted from each other in the PDF-Export.

Symbols in PDF-Export File:

Symbols in QGIS


texts: How to work with raster cell sizes when using Focal Statistics?

texts:
I'm trying to calculate the percent area of a certain land cover type within a 5000m window but my results aren't as expected. I've never used focal statistics before but I'm almost certain the results retrieved are not correct and it has something to do with the raster cell size.
I downloaded a landcover raster from Google Earth Engine and the cell size is 0.000808483755707569 for both X and Y.
In the focal statistics tool I used the following parameters:
Neighborhood: Rectangle
Width: 5000
Height: 5000
Unit Type: Cell
I retrieved the following output:

Is this a correct output?


texts: map.on('click') event disables the double click zoom

texts:
I have a map that has an onclick event fired and displays a marker where the user has clicked. Problem is that this event seems to overrun the double click event. Because of this, when I double click the map, the map puts a marker and doesn't zoom. Anyone has an idea how I can solve this?

openmap = L.mapbox.map('map', 'julievelghe.nmjoj53p',{doubleClickZoom: true}).setView([49.152969, 6.473007], 5);

//on map click, create a marker
openmap.on('click',function(event){
    var coordinates = event.latlng;
    placeMarker(coordinates);
    alert("Single click");
});


texts:
Problems trying to reproject the DEM raster layer onto the city layer:
The base layer (CITY) is using NAD_1983 StatePlane Georgia West FIPS1002 feet (WKID:2240) projected CRS with a geographic CS: GCS_NAD1983 (WKID:4269)
Raster Dataset extent: Elevation surface for Cobb County, GA including City of Austell.
West: -89.857908 East:-89.718310
North 40.249497 South: 40.090210
What I have tried :

I tried raster functions -> reproject
I tried creating an unknown_CS layer, overwrite the original DEM coordinates then define it again with the town coordinates using DEFINE PROJECTION
I tried "project raster" tool

Yet the layers still do not align properly, though when I view the meta of the DEM file, the stated projected CRS and geographic CS are the same. What is causing this issue?


texts: I have a problem in Qgis where all my labels are doubled on every single item (see image). I have checked the attributes list and everything seems normal. How do I show the labels properly?


texts:
I am using ESRI Sliverlight API 3.1 using layer definition in dynamic service layer, when I create a definition ID in ('1','2','3'......'1000001') . This definition is not getting applied, when I pass big string in REST unable to perform the query. But I need to set a filter based on ID, but I cannot number of ids. How to resolve the issue?

                LayerDefinition def = new LayerDefinition();
                def.LayerID = layerId;
                def.Definition = "ID" + " in " + idString;
                dynamicServiceLayer.LayerDefinitions.Add(def);


texts: Geocoding in ARCGIS desktop with online basemaps

texts:
I would like to have your guidance about possible usage of online base maps in ArcGIS.

In ArcGIS Desktop there is a possibility to have access to some online maps (Satelite maps, Street Msps,...) provided by ESRI as backgroud info.

I would like to know if there is a possibility to do geocoding using this maps for some points.


texts:
I am trying to visualize live data that is updated every 5 seconds coming from a GeoJSON layer served by GeoServer.

What is the best way to do this on Cesium?

I saw that the GeoJsonDatasources in Cesium is intended only for static data, does this mean I have to convert my GeoJSON to CZML to be able to use the CzmlDataSource or it's a better choice to build a custom datasource and update the entities using the GeoJSON layer that is being served by GeoServer?


texts: Heads up before trying Mapbox vector tiles services flexibility with OpenLayers 3

texts:
if I load Mapbox vector tiles services in my map using OpenLayers 3, can I perform the functions like on click popup or on hover popup, highlight etc?

Is it advisable and mature enough? As I saw this example it has a warning of browser becoming unresponsive. Has anyone tried to load say ten vector tiles in one map?


texts: How to reshape an array of points into an image

texts:
I have a numpy array of shape (1964, 3), in the following format:

[
   [ 0.3945, 52.700, 0.8 ],
   [ 0.3788, 52.121, 0.3 ]
]
etc.

Each nested array represents a point in the format of [longitude, latitude, value]

I am trying to do a points.reshape() in order to convert this array into an image however am unable to calculate the columns and rows for the image.

It might be that I'm missing something incredibly simple here but can anyone suggest a solution to this issue?

I also have a pickle file of the array I'm working with so I'm happy to send if anyone wants.


texts: Adding Permalink attribution to the leaflet map in order to reach another URL permalink

texts:
Continuing the issue opened here:

Link to another URL permalink in Leaflet

I found some code, that may help me to establish a link to another url permalink from the leaflet map:

https://raw.githubusercontent.com/MapBBCode/mapbbcode/master/src/controls/PermalinkAttribution.js

The library is the "PermalinkAttribution".
Inside this file important is:

 for (var i in this._attributions) {
        if (this._attributions[i]) {
            // make permalink for openstreetmap
            if( i.indexOf('/openstreetmap.org') > 0 ||
 i.indexOf('/www.openstreetmap.org') > 0 ) {
                var permalink = 'http://www.openstreetmap.org/#map=
 {zoom}/{lat}/{lon}';
                i = i.replace(/(['"])http[^'"]+openstreetmap.org[^'"]*
  (['"])/, '$1' + permalink + '$2');
                if( this._map.options.attributionEditLink ) {
                    var editlink = permalink.replace('#', 'edit#');
                    i = i.replace(/(openstreetmap.org[^'"]*(['"])[^>]*>[^<]+
  <\/a>)/, '$1 (Edit)');
                }
            }
            var latlng = this._map.getCenter();
            i = i.replace(/\{zoom\}/g, this._map.getZoom()).replace(/\
  {lat\}/g, L.Util.formatNum(latlng.lat, 4)).replace(/\{lon\}/g,
  L.Util.formatNum(latlng.lng, 4));
            attribs.push(i);
        }
    }

but so far major problem for me was how to add it on the map.

I tried:

   L.Control.Attribution.extend({}).addTo(map);

   L.Control.PermalinkAttribution{}).addTo(map);

but nothing works.

Console says:

throw new Error('Map container not found.');

and I have no button on  my map shown.

Does anyone knows how to attach the permalink attribution to the leaflet map?


texts:
Is there a way to serve mapbox GL  gzipped geojson files?
I have style.json file that pulls multiple geojson files and it is pretty heavy.


texts:
I have satellite image files with metadata generated from ArcGIS:

Raster Download Information

Coordinate system ID WGS84:  4326
Top edge WGS84:  35.0003681085064 Degree
Bottom edge WGS84:  34.991149768818 Degree
Left edge WGS84:  -80.4943741887217 Degree
Right edge WGS84:  -80.4841387811576 Degree

I want to map this back to Google map's lat and longtitude, so that I can have overlay of street on the image. I used topedge for lat and left edge for long, but get a different place. How can I find the transformation?


texts: Using Python add-in wizard with ArcGIS 10.3 for Desktop?

texts:
I'd like to create an add-in for ArcMap (ArcGIS 10.3) using Python. Unfortunately, the wizard does not seem to have been updated for 10.3, which is the only version Esri provides for download at present (AFAIK). I don't have any other version installed, and never did because I'm starting fresh on a new PC.

Is there a wizard for 10.3?

If not, is there good documentation about what the wizard did for 10.1 and 10.2? If it only creates an XML document supplying the metadata and the name of the Python class to be instantiated on button-push, I'm sure I can replicate that. I just can't find any detailed information from ESRI; the help documents discuss the wizard as though it is the only option.


texts: How to position direction arrows on lines' sides with Geoserver?

texts:
I'm using Geoserver 2.6.2. I need to represent direction arrows for line features. I have managed to position such marks on lines centroid but I need to position these arrows on the side of the lines. The perpendicularOffset option used for line labels would be the solution I'm looking for but it is restricted to labels styling. How would you code your SLD in order to achieve that?

Below is the code I'm using so far to position arrows on lines:

    <Rule>
      <Name>Lines</Name>
      <Title>Lines</Title>
      <LineSymbolizer>
        <Stroke>
          <CssParameter name="stroke">#000000</CssParameter>
          <CssParameter name="stroke-width">1</CssParameter>
        </Stroke>
      </LineSymbolizer>
      <PointSymbolizer>
        <Geometry>
          <ogc:Function name="centroid">
            <ogc:PropertyName>the_geom</ogc:PropertyName>
          </ogc:Function>
        </Geometry>
        <Graphic>
          <Mark>
            <WellKnownName>shape://oarrow</WellKnownName>
            <Fill>
              <CssParameter name="fill-opacity">0</CssParameter>
            </Fill>
            <Stroke>
              <CssParameter name="stroke-width">1</CssParameter>
            </Stroke>
          </Mark>
          <Size>20</Size>
          <Rotation>
            <ogc:Function name="endAngle">
              <ogc:PropertyName>the_geom</ogc:PropertyName>
            </ogc:Function>
          </Rotation>
        </Graphic>
      </PointSymbolizer>
    </Rule>


texts: Counting trees in a vector shapefile layer using PyQGIS

texts:
[The result (shp file) My goal is to count trees in shapefile and determine whether it is olive tree or not.
For now I have .shp files and I can't count how many trees in it.
Also is it possible to export every tree to a shp file by itself?
This is my code:
from osgeo import gdal, ogr
import sys
import os
import gl   ob
from osgeo import ogr
from osgeo import gdal
Select_folder = "path"

x=1
while x < 6:

        fn = Select_folder + "/" + str(x) + ".tif"
        fi = QFileInfo(fn)
        fileName = fi.baseName()
        gdal.UseExceptions()
        os.chdir("path")
        fileName = Select_folder + "/" + str(x) + ".tif"
        src_ds = gdal.Open(fileName)
        if src_ds is None:
            print('Unable to open %s' % src_fileName)
            sys.exit(1)
        srcband = src_ds.GetRasterBand(1)
        dst_layername = "PolyFtr" + str(x)
        drv = ogr.GetDriverByName("ESRI Shapefile")
        dst_ds = drv.CreateDataSource(dst_layername + ".shp")
        dst_layer = dst_ds.CreateLayer(dst_layername, srs = None)
        newField = ogr.FieldDefn('Area', ogr.OFTInteger)
        dst_layer.CreateField(newField)
        gdal.Polygonize(srcband, None, dst_layer, 0, [],callback=None )
        x+=1

y=1
while y < 6:
        #rlayer = iface.addRasterLayer(fn,fileName)
        yn = Select_folder + "/PolyFtr" + str(y) + ".shp"
        yi = QFileInfo(yn)
        yileName = yi.baseName()
        layer = iface.addVectorLayer(yn, yileName , "ogr")
        y+=1


texts:
How can the mid-point of a polyline be calculated using the geoprocessing framework, e.g. in a Python script?

The polyline.centroid property returns The true centroid if it is within or on the feature; otherwise, the label point is returned. The centroid is rarely located on non-straight lines, which is useless for my purposes.

The Feature Vertices to Points tool has a midpoint option but this requires ArcInfo, which I don't currently have.

An option could be to add measures to the polylines and create a route event 50% along the line.

Another workaround is to use the Calculate Geometry option in ArcMap, but ideally I need to automate this process in a script.

Any better/faster suggestions?

I am limited to ArcGIS 10.0 (without access to an ArcInfo level license) for the moment.


texts: Creating Topographic Position Index in SAGA GIS?

texts:
I am trying to create a TPI from a DEM tiff file in SAGA but keep getting the following result:

I'm just not sure what I'm doing wrong.


texts:
I have a layer with a column conaining population values. I want to write an expression to compute in a new field the percentage of population each polygon has in relation to the total population.

I know I can use basic stats to get the SUM of area column, but I want an automatic process to repeat over several fields).


texts: Is it possible to put the projection of a map of a layout in a text label in a layout in QGIS?

texts:
I would like to be able to have a text label in a print designer with the information of the projection of a map in QGIS. That way, if I change the projection of a project, I would change that information within the layout.

I think I've seen a variable map_src, but I can not find a way to use it.


texts: Keeping properties secure within mapbox-gl

texts:
Simply put, how could one develop a mapbox-gl application where the properties are secured ?

For example, take a look at this simple mapbox-gl map with properties:

http://jsbin.com/hovawa/edit?html,js,console,output

In the code I have property called "secretproperty" with value "abc123".

Is it possible to somehow hide/encode this if someone attempted to read the property info in the browser's debug mode?
What would be the best approach?


texts: Arcgis table and Excel sheet dynamic linking

texts:
How i can link an ArcGIS table with Excel sheet, so the table updated automatically when i update the excel sheet?

Also, can i have sub-table in ArcGIS?


texts: Geoserver slow rendering un-cached layers

texts:
I have GeoServer connected to a SQL Server datastore.  I have a map with 5 Layers, the layers have about 200,000 features total.  The features are styled as simple points or lines - not many lines though. Due to the nature of the data - it changes alot and needs to be instantly current - I don't have caching on.

Geoserver is on a windows server, running on Tomcat. This is all being rendered in the browser using OpenLayers.

I have added some parameters to the CATALINA_OPTS environment variable -

CATALINA_OPTS=-Xms128m -Xmx2048M -XX:SoftRefLRUPolicyMSPerMB=36000

and installed the two recommended Geoserver and JAI and ImageIO extensions.

I have tried rendering as TileWMS and ImageWMS.

Peformance is sub-optimal at best.  It really looks like SQL server returns data very fast, and the bottleneck is GeoServer but I have more questions than answers at this point.  It really doesn't look like the Tomcat service is getting used hard watching it's process in Task Manager.

How can I optimally set up Geoserver to give the best performance? Is there a best practice guide (besides http://docs.geoserver.org/latest/en/user/production/index.html) to point me in the right direction?


texts:
I'm trying to add a CartoDB basemap to my OpenLayers map but I can't get it to work. Here is the code I use to declare the basemap variable:
const carto = new ol.layer.Tile({
    source: new ol.source.XYZ({
        url:'http://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png',
    })
});

Any idea where the problem comes from?


texts: My UpdateCursor is not working. The attribute table of my new shapefile is not updating

texts:
I am new to arcpy. I am trying to do a few things with a shapefile. First I created a new shapefile with specific dataset that I wanted to work on. So far everything works until the UpdateCursor. When it runs, it doesn't update the attribute table of the output file (outfc). Am I missing something?

Here is what I'm trying to do:

All of the values that are in the “SPCOMP” field should be put into related individual species fields, and leave the rest of species fields to 0. For example, if the “SPCOMP” = “Pj 60Pt 30Bw 10”, that means that Pj = 60, Pt = 30, Bw = 10. Then 60 should go under “PJ”, 30 goes under “PT”, and 10 goes under “BW”.

Here is the code:

import arcpy
from arcpy import env

env.overwriteOutput = True

env.workspace = 'G:\\School\\GGR386\\FRMG_Coding_Exercise'

infc = 'ftg_true_all.shp'

outfc = 'G:\\School\\GGR386\\FRMG_Coding_Exercise\\New\\ftg_true_all_change.shp'

where_clause = """ "TRTMTHD1" = 'NATURAL' OR "TRTMTHD1" = 'PLANT' """

arcpy.Select_analysis(infc, outfc, where_clause)

in_field = 'SPCOMP'
fields = arcpy.ListFields(outfc)
field_list = []

sc = arcpy.SearchCursor(outfc, [in_field])

for row in sc:
    s = str(row.getValue(in_field))
    s_stripped_cap = s.replace(' ', '').upper()
    i, j = 0, 0
    spcomp = []

    if s_stripped_cap[0].isdigit() == True:
        data_type = 'int'
    else:
        data_type = 'alphabet'
    for char in s_stripped_cap:
        if char.isdigit() == True and data_type != 'int':
            spcomp.append(s_stripped_cap[i:j])
            i = j
            data_type = 'int'
        elif char.isalpha() == True and data_type != 'alphabet':
            spcomp.append(s_stripped_cap[i:j])
            i = j
            data_type = 'alphabet'
        j += 1

    for f_name in fields:
        if f_name.type != "String" and len(f_name.name) <= 2:
            field_list.append(f_name.name)

        spuc = arcpy.UpdateCursor(outfc, [field_list])
        for r in spuc:
            for i in range(len(spcomp)):
                if spcomp[i].isalpha() and spcomp[i] == r:
                    sp_num = int(spcomp[i + 1])
                    r[0] = sp_num
                    spuc.updateRow(r)

del row
del sc
del r
del spuc


texts: What are useful websites for scouting GIS related jobs?


texts:
I want to create a standalone pyQGIS script in python editor but I have importing error on Processing

import sys
from qgis.core import *
sys.path.append('/usr/share/qgis/python/plugins')
import Processing

I turn for similar question and all suggestion import path for Processing plugins and I add this path and I have error again.
How to fix that on Ubuntu 16.04


texts: I have upgraded from 10.0 to ArcGIS 10.1.  I am trying to open an 10.0 mxd file in 10.1 and it appears to have lost any additions I had made to the table plus its' x and y coordinates.  Can you successfully migrate ArcGIS 10.0 mxd to 10.1?


texts:
I need update/refresh/repaint the layer on the map when i click on the button refresh.
I remove the layer and add it again but it not working and points on the map are not displaying.
Is there an other way to refresh data.

private void showMap(File rasterFile, final File shpFile) throws Exception{

    AbstractGridFormat format = GridFormatFinder.findFormat( rasterFile );
    reader = format.getReader(rasterFile);

    final Style rasterStyle = createRGBStyle();

    // Connect to the shapefile
    FileDataStore dataStore = FileDataStoreFinder.getDataStore(shpFile);
    SimpleFeatureSource shapefileSource = dataStore.getFeatureSource();

    Style shpStyle = SLD.createSimpleStyle(dataStore.getFeatureSource().getSchema(), Color.RED);

    // Set up a MapContent with the two layers
    final MapContent map = new MapContent();
    map.setTitle("The Map");

    Layer rasterLayer = new GridReaderLayer(reader, rasterStyle);
    map.addLayer(rasterLayer);

    final Layer shpLayer = new FeatureLayer(shapefileSource, shpStyle);
    map.addLayer(shpLayer);

    // Create the mappan.
    final JMapPane mapPane = new JMapPane(map);

    jInternalFrame1.setLayout(new BorderLayout());
    jInternalFrame1.add(mapPane, BorderLayout.CENTER);

    JMenuBar menuBar = new JMenuBar();
    jInternalFrame1.setJMenuBar(menuBar);
    JMenu menu = new JMenu("Map");
    menuBar.add(menu);

    menu.add( new SafeAction("Refresh") {
        public void action(ActionEvent e) throws Throwable {
            map.removeLayer(shpLayer);
            // Connect to the shapefile
            FileDataStore dataStore = FileDataStoreFinder.getDataStore(shpFile);
            SimpleFeatureSource shapefileSource = dataStore
            .getFeatureSource();

            Layer shpLayerNew = new FeatureLayer(shapefileSource, shpStyle);

            map.addLayer(shpLayerNew);
        }
    });
}


texts: Changing default DPI for Export Map to Image in QGIS?

texts:
In QGIS 3, is there a way, perhaps a configuration file, to change the default DPI for doing Export Map to Image from the Project \Import/Export menu?

I use this function a lot, but the default DPI, at 96 dpi, is too low, and every time I need to manually change it to 300 dpi.

This is with QGIS 3.4, Ubuntu 18.04 LTS.


texts: Simpler Method to Map Fields in Spatial Join

texts:
I've cobbled together a script, part of which creates field mapping to sum fields in a spatial join.  It took 13 lines of codes just to define which fields I want summed in the output.  Is there a simpler/more elegant way to define and sum output fields?

import arcpy, os

ws = arcpy.GetParameterAsText(0)
ws2 = arcpy.GetParameterAsText(1)
outws = arcpy.GetParameterAsText(2)

arcpy.env.workspace = ws

# List all feature classes in ws
fcs = arcpy.ListFeatureClasses()

for fc in fcs:
    # Create full paths for pair1 and pair2
    pair1 = os.path.join(ws, fc)
    pair2 = os.path.join(ws2, fc)
    outfc = os.path.join(outws, "Joined_" + fc)

    #Create Field Mappings
    fieldMappings = arcpy.FieldMappings() #Create field mappings object

    fieldNamesToSum = ['Min_MA_NEW', 'Min_01_NEW', 'Min_02_NEW', 'Min_03_NEW', 'Min_04_NEW', 'Min_05_NEW', 'Min_06_NEW', 'Min_07_NEW', 'Min_08_NEW', 'Min_09_New', 'Min_10_NEW', 'Min_11_NEW', 'Min_12_NEW'] #Field names to sum

    # This creates a field mapping that only outputs the sums of the joined fields in fieldNamesToSum - the range number must equal the number of fields in fieldNamesToSum
    for i in range(13):
        fm_type = arcpy.FieldMap()
        fm_type.addInputField(pair2, fieldNamesToSum[i])
        outField = fm_type.outputField
        outField.name = fieldNamesToSum[i]
        fm_type.outputField = outField
        fieldMappings.addFieldMap(fm_type) #Sequentially add 13 FieldObjects
        fieldIndex = fieldMappings.findFieldMapIndex(fieldNamesToSum[i])
        fieldMap = fieldMappings.getFieldMap(fieldIndex)
        fieldMap.mergeRule = 'Sum'
        fieldMappings.replaceFieldMap(fieldIndex, fieldMap)

    arcpy.SpatialJoin_analysis (pair1, pair2, outfc, "JOIN_ONE_TO_ONE", "KEEP_ALL", fieldMappings, "SHARE_A_LINE_SEGMENT_WITH")


texts: Clipping line with overlaying line in QGIS

texts:
I have two layers of line which overlap as seen in the picture.
I would like to cut out the first picture with the other as an overlap and then merge them together to get one single layer.
Or is there another way to do this easier?


texts: QGIS programmatically change outline colour of circles

texts:
I have a vector layer with a "Simple Marker" for the style. I see I can create Data Defined properties for both the colour and the outline colours.

In Python, I'm able to change the marker's base colour (for the data defined property), but I can't seem to find how to change the outline's colour.

Here's what I have:

symbols = layer.rendererV2().symbols()
symbol = symbols[0]

expr = 'color_rgba(255,0,0,0)'
symbol.symbolLayer(0).setDataDefinedProperty('color', expr)

# I would expect this to work, but it doesn't
symbol.symbolLayer(0).setDataDefinedProperty('outlineColor', expr)

Any ideas?


texts: Do not show a point symbol for a point if closer than nn meters to another point in same layer

texts:
I have a layer consisting of points (big stones), and there are too many points to be able to show them all in the scale I am using in my print composer.

What is the most efficient, and preferably least processor heavy, way to tell QGIS to only show a point, if there is no other point in the same layer which is closer than, say, 50 meters, from another point in that layer. There is in this case no preference which stone to show and which to not show, ie I have no data telling for instance which stones are bigger or more important.


texts: Choosing a platform for developing opensource GIS web application.

texts:
I am from a non programming background. I started learning GIS opensource web applications from last 3 moths. Now I am using only notepad++ for writing codes and I am familiar with javascript, HTML and litle bit PHP.
My question is,
              1.Can I develop a application without java or asp.net.
              2.If yes , text editor (notepad++) is sufficient or I have to go with eclips, netbeans or visualstudio.


texts: How to clip an Openlayers OSM Background to a Vector Boundary?

texts: I'm new to QGIS. In my project i have two layers: one is my vector (city boundaries) and the other is the OSM (obtained with the plugin). How can I cut the OSM along the edges of the vector? I want to obtain only a region of the map, to save it as an image.


texts:
I've encountered some problems running my old 1.8 project in 2.0. I would like to reinstall QGIS 1.8 but seems that there are no more 1.8 installer  in KingChaos. In QGIS "all releases" there isn't.

Where I can find the old release for Mac user?


texts: Where does QGIS write its error log messages?

texts:
I want to read and copypaste the error log for the crashes I am experiencing in QGIS 2.0. The error I am getting freezes QGIS so I have to ctrl-alt-del. I can see the last message in the panel, but can't copy it or scroll back to see the earlier messages. When I re-open it, the log message panel has reset to zero (QGIS starting, etc), so all the transform error messages I want are gone.

Are the older log messages written to disk somewhere?


texts: Packing polygons without intersection and empty places

texts:
How can I pack polygons without empty places and intersections?
Please tell me everything that might help me (API, algorithms, etc.)

I'm using:

USPS API to get carrier routes.
dbscan.js to generate clusters from points.
hull.js to generate polygons.
leaflet to show map

Example (I have):

Example (I need):


texts: Speeding up Arcpy script with lots of for loops and cursors?

texts:
I'm running this chunk of code (below) inside a larger script. This section is supposed to find the time for each boat (mmsi) that is closest before and closest after a single detection time. It first finds all the boat times before the detection time and puts them into a fc, then it does the same with all the boat times after the detection time. It then finds the max and min (respectively) of those two fcs and put them into yet another fc. However, this chunk of code is a) messy and b) taking forever to run.

I know from reading about speeding up python scripts that for-loops and  cursors can slow things down. However, because of the way they are all nested in code, I can't figure out how to rewrite this in a simpler, faster way. Any thoughts?

arcpy.env.workspace = r'L:\gathr\indonesia\Sara\Date20160104\Detections_UTM.gdb'
workspace = r'L:\gathr\indonesia\Sara\Date20160104\Detections_UTM.gdb'
Geodatabase = arcpy.ListFeatureClasses()
detection_datalist = []
print Geodatabase
for fc in Geodatabase:
    fc_name = str(fc)
    print fc_name
    detection_name = fc_name
    fc_path = os.path.join(workspace, fc_name)
    cursor = arcpy.SearchCursor(fc_path)
    for row in cursor:
        time1 = row.getValue("Date_Mscan")
        time2 = time1.strftime("%Y-%m-%d %H:%M:%S")
        print time2
        time = datetime.strptime(time2,"%Y-%m-%d %H:%M:%S")
        print time
    del cursor
    arcpy.env.workspace = r'L:\gathr\indonesia\Sara\Date20160104\fc_processed.gdb'
    Geodatabase = arcpy.ListFeatureClasses()
    for fc in Geodatabase:
        fc_str = str(fc)
        suffix = detection_name
        if fc_str.endswith(suffix):
            where_clause = '"ts_pos_utc_Converted" <= date' + "'%s'" %time
            before = os.path.join(r'L:\gathr\indonesia\Sara\Date20160104\fc_processed.gdb', fc + "_before")
            arcpy.Select_analysis(fc, before, where_clause)

            listdates = []
            rows = arcpy.SearchCursor(before)
            row = rows.next()
            while row:
                value = row.ts_pos_utc_Converted
                print value
                listdates.append(value)
                row = rows.next()
            print listdates
            try:
                maxdate = max(listdates)
            except:
                maxdate = None
            print maxdate

            after = os.path.join(r'L:\gathr\indonesia\Sara\Date20160104\fc_processed.gdb', fc + "_after")
            where_clause = '"ts_pos_utc_Converted" > date' + "'%s'" %time
            arcpy.Select_analysis(fc, after, where_clause)

            listdates = []
            rows = arcpy.SearchCursor(after)
            row = rows.next()
            while row:
                value = row.ts_pos_utc_Converted
                print value
                listdates.append(value)
                row = rows.next()

            print listdates
            try:
                mindate = min(listdates)
            except:
                mindate = None
            print mindate
            new_fc = os.path.join(r'L:\gathr\indonesia\Sara\Date20160104\fc_processed.gdb', fc + "_final")
            layer = os.path.join(r'L:\gathr\indonesia\Sara\Date20160104\fc_processed.gdb', fc)
            arcpy.CopyFeatures_management(layer, new_fc)
            with arcpy.da.UpdateCursor(new_fc,"ts_pos_utc_Converted") as cursor:
                for row in cursor:
                    if row[0] not in [mindate, maxdate]:
                        cursor.deleteRow()

del row, cursor


texts: Sample size for variogram modelling using the geoR package in R

texts:
I am trying to use the geoR package in R to model the variogram for a dataset of 15 000 sample points of bathymetry data. I intend to do a random split to create test (20%) and training samples (80% of the samples), so the accuracy of my final predictions can be assessed.

However, when using all the sample points (15000) I get the following error:

variog: computing omnidirectional variogram
Error: cannot allocate vector of size 427.9 Mb

The error is associated with the large number of samples, because when I reduce the sample size to about 5000, it runs perfectly. I have tried using a machine with 32G of RAM on the 15 000 samples, to no avail.

Would it be practical to use a sample of the data solely for the purposes of the semi-variogram modelling? If not, what would be the best approach to dealing with this problem? Any suggestions would be welcomed.


texts: Arcpy Cell Statistics - maximum number of input rasters

texts:
I am running cell statistics in arcpy. Essentially I am summing together a folder of grd files to create a new grd representing their sum

This will eventually take place on a folder of 20,000 raster grds.

Is there a maximum number of raster datasets that can be processed by the cellStatistics tool?


texts: Making ee.data.authenticationViaAuth always success

texts:
I am new to Google earth engine. I am using earth engine javascript api (ee_api_js.js).

var onImmediateFailed = function() {
    ee.data.authenticateViaPopup(function() {
      runAnalysis();
    });
};
ee.data.authenticateViaOauth(CLIENT_ID, runAnalysis, null, null, onImmediateFailed);

When I run it the first time on any browser it always goes to "onImmediateFailed". I am not sure what settings needs to be done in order to always succeed while executing this.

I am creating a web app where users will navigate to earth engine functions and I don't want user to have to login to their Google account.


texts: Calling PyQGIS functions from external editors (Linux)

texts:
Is it possible to integrate external Python editors (such as KDevelop) with QGIS, so that it would be possible to run functions in qgis.core, qgis.utils etc. outside of the QGIS Python Console?

Following the guidelines on the QGIS website (http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/intro.html), I tried this, but it just returns 1 and nothing else:

import sys
sys.path.append('/usr/share/qgis/python')
import qgis.core
import qgis.utils

app = qgis.core.QgsApplication([], True)
qgis.core.QgsApplication.initQgis()
qgis.utils.iface.addVectorLayer("testing.shp", "anewlayer", "ogr")
aLayer = qgis.utils.iface.activeLayer()
print aLayer.name()

Such as:

$ LD_LIBRARY_PATH=/usr/lib64/qgis/ python qgis-test.py && echo "OK" || echo "Died"
Died

I am running openSUSE Tumbleweed, 64-bit.


texts: KML to Shapefile using GDAL (LIBKML) - MultiGeometry data not transforming

texts:
I am trying to convert a KML file to Shapefile using GDAL (ogr2ogr).
Regular KML files convert fine but KML files containing MultiGeometry features (eg: Geometry containing both Point and Polygon) do not get transformed.

The output is simply a blank shapefile containing no shapes or attributes.

I have tried this command:

ogr2ogr -overwrite -f "ESRI Shapefile" -where OGR_GEOMETRY='MultiGeometry' <ouputfolder> <input kml>

For an example, KML file which has MultiGeometry, you can see the us_states.kml hosted by Google.

US States KML - Google.com

Other info: I am using 'GDAL 1.10.1, released 2013/08/26' for windows.


texts:
I have a list of data points with timestamps in the format:

row: Object ID, timestamp

I would like to find the minimum timestamp grouped for each object id, something like:
minimum(timestamp, object ID)

this format works when I have other data formats. Is this because my timestamp is in the format Date&Time?
Any suggestions?


texts:
I have downloaded a historical USGS map that is suppose to be a geoPDF. My understanding is that QGIS cannot directly read geoPDF's. I have not had any luck in converting it to a geoTIF in order to bring into QGIS.

Has anyone had any success in doing this? Is there another way that
someone knows to use the geoPDF with QGIS.
I converted the geoPDF to a TIF and attempted to use the
georeferencer plug-in. Everytime I try doing it te results are not
correct.


texts:
Need ideas.

I have a  layer of city shapefile with zones:

I have a layer with FCD data:

My goal/question is:

Is it possible to define a model where I will be able to choose any 2 zone and get points started and ended (not just passed through(!)) between these 2 given zones?

For example: id_easy - d1ed640cca started and finished trip between given zones.

FCD data format (very small sample):

id_easy    ordinal  latitude    longitude   epoch             day_of_week   geometry
d1ed640cca       1  11.0589     3.6907      04-06-16 00:48    Saturday      POINT (3.6907 11.0589)
745b0999ce       1  11.0591     3.6911      04-06-16 00:40    Saturday      POINT (7.6911 11.0591)


texts: I am working on Land Use/Land Cover Classification of the Hindon river(India) watershed extracted using Aster DEM. Landsat 8 scenes had been used and mosaicked for the required AOI using the Watershed shapefile. I am facing a lot of troubles in the Unsupervised Classification, ISODATA method as I am getting numerous mixed pixels which are difficult to label. I had also tried histogram equalization before clustering but the recoding & labelling issues persist. Kindly suggest me some procedure to overcome the difficulties faced in labelling the classes.


texts: How to get the co-ordinates of the clicked point using 'ol/Draw' in OpenLayers 6?

texts:
I am trying to create a measure tool using OpenLayers 6. I did get examples to measure distance and area using the ol/Interaction/Draw and ol/geom library classes. I want also to measure the coordinates of the clicked point. But I am unable to use the draw class when using "type: 'Point'"

I have come across examples in which this is implemented using the 'singleclick' event trigger on the specified map. But I am looking forward to do the same task using the latter mentioned method.
How can I do this?

I have implemented measure distance and area using the below code

mapHandlers() {
    const that = this;
    let listener;
    this.draw.on('drawstart', (evt) => {
      // set sketch
      that.sketch = evt.feature;
      let tooltipCoord = evt.coordinate;
      listener = that.sketch.getGeometry().on('change', (evt) => {
        const geom = evt.target;
        let output;
        if (geom instanceof Polygon) {
          output = that.formatArea(geom);
          this.measureArea = (getArea(geom, {projection: 'EPSG:4326'}));
          tooltipCoord = geom.getInteriorPoint().getCoordinates();
        } else if (geom instanceof LineString) {
          output = that.formatLength(geom);
          this.measureDistance = (getLength(geom, {projection: 'EPSG:4326'}));
          tooltipCoord = geom.getLastCoordinate();
        }
        that.measureTooltipElement.innerHTML = output;
        that.measureTooltip.setPosition(tooltipCoord);
      });
    });

    this.draw.on('drawend', () => {
      that.measureTooltipElement.id = 'ol-tool';
      that.measureTooltipElement.className = 'ol-tooltip ol-tooltip-static';
      that.measureTooltip.setOffset([0, -7]);
      // unset sketch
      that.sketch = null;
      // unset tooltip so that a new one can be created
      that.measureTooltipElement = null;
      that.createMeasureTooltip();
      unByKey(listener);
    });
  }

I am referring to, measure example in openlayers


texts:
Is it possible for QGIS to read ESRI ST_GEOMETRY in an Oracle database? The ESRI geometry type looks like this:

SDE.ST_GEOMETRY(260,61,491455.9942,5456922.8891,492388.3189,5457507.9101,NULL,NULL,NULL,NULL,0,863.557608125545,26910,(BLOB))

Alternatively, is there an SDE function or group of functions that would represent the ST_GEOMETRY as SDO_GEOMETRY? Then I could make a view for the ST_GEOMETRY and view it in QGIS and other non-ESRI software.


texts: Leaflet polygon (multi polygon) to CQL equivelent

texts:
I am using leaflet with geoserver. I want to set up a CQL filter to do an intersect but having an issue with getting the geometry from leaflet polygon to cql polygon.

For instance, the feature has the following.

gemoetry.coordinates is something like this when output as a string.

-90.89,46.67,-90.89,46.67,-90.89,46.67,-90.89,46.67,-90.89,46.67

CQL would need it in this format.

POLYGON((-90.89 46.67,-90.89 46.67,-90.89 46.67,-90.89 46.67,-90.89 46.67))

This would then be used to do something like this.

cql_filter=INTERSECTS(Geometry, POLYGON((-90.89 46.67,-90.89 46.67,-90.89 46.67,-90.89 46.67,-90.89 46.67))

However, it doesn't seem like leaflet has a toString or toCQL to convert a leaflet polygon to a CQL polygon. Now, one thing that comes to mind is writing a custom parser to go through and generate the CQL polygon needed to do the filter. This doesn't seem ideal to me and not sure how it would handle multipolygons. Is there an easy way to do this is leaflet?

UPDATE: Based on the answer, I used the wicket with the following code.

var wkt = new Wkt.Wkt();
wkt.fromObject(object.geometry)
var filter = "INTERSECTS(Geometry," + wkt.write() + ")";


texts: comparing and finding inaccuracies in two raster layers

texts:
I have two raster layers with urban cells and rural cells. I used an overlay function to compare these two raster layers and see the differences. Before overlaying I gave urban cell  the value 1 and rural cell the value 0. After overlaying I get the result: 0,1,2,NA:
Clearly, 2 and 0 mean that the urban and rural cells match and 1 mean it doesn't. however just looking at 1 i cannot see if the error was occured due to urban cell or rural cell.
Does anyone know how can I see where the inaccuracy in this raster layer is if i take one of the raster layers as reference?

b
class       : RasterLayer
dimensions  : 190, 333, 63270  (nrow, ncol, ncell)
resolution  : 0.008344257, 0.008344257  (x, y)
extent      : 12.10615, 14.88479, 51.68836, 53.27377  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

 c
class       : RasterLayer
dimensions  : 190, 333, 63270  (nrow, ncol, ncell)
resolution  : 0.008344257, 0.008344257  (x, y)
extent      : 12.10615, 14.88479, 51.68836, 53.27377  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

 uc <- overlay(c, b, fun=function(x,y){return(x+y)})
 urx <- overlay(c, b, fun=function(x,y){return(x*y)})

 cellStats(uc, stat = "sum",na.rm = TRUE)
 unique(getValues(uc))
 [1] NA  0  1  2


texts: How to access existing graphs in arcGIS 10

texts:
I am trying to access an existing graph in arcGIS version 10.  The student workbook I am using was written for version 9.3.  Version 10 does not have "tools" option in the toolbar. The instructions guide students to the 'Tools' option.
I will be teaching students to use arcGIS using 'Mapping our World' student workbooks (written for version 9.3) using version 10.


texts: Inverse distance weighted interpolation for categorical and/or binary datasets

texts:
I want to perform an Inverse distance weighted (IDW) interpolation to a series of points. However my data set can be categorical or binary and just like the IDW I want to interpolate the new values (categorical or binary) based on the highest frequency of closest points values surrounding the prediction location:

A if the majority of the closest points are A
C if the majority of the closest points are C
...

In addition I want to avoid using close points that are behind a barrier, that can be a wall or just an elevation.
This is how the problem/dataset looks like:

The grey points is the point I want to interpolate and the points behind the barriers can't be used in the process. The output of the interpolation process would be A or B, in this case probably A, since the closest points are all A.
What I am looking for is some kind of work-flow that would allow me to do this in ArcGIS, if possible I would prefer R (if anyone has any experience in doing this kind of processes in R). This is just the example how my problem looks like, I still have to implement it to a bigger extent and to many layers of points. Bear in mind that I have little experience with interpolation.


texts:
I can edit the line vertices I create in Collector, but can't figure out how to edit the vertices directly in ArcGIS online.

I can change the length of the line, but not the vertices in ArcGIS online


texts: Field not updating correctly by updatecursor() method

texts:
in my work, I have to find the features having duplicate values of field "ID" and leave just one identical value by updating blank values in rest of the features.

My code for an add-in is:

import arcpy
from arcpy import env
import pythonaddins

class ButtonClass1(object):
    """Implementation for Remove_Duplicate_ID_addin.button (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        mxd = arcpy.mapping.MapDocument("CURRENT")
        layers = arcpy.mapping.ListLayers(mxd)
        for layer in layers:
            desc = arcpy.Describe(layer)
            path1 = str(desc.path) + r"\outstat"

            env.workspace = desc.path
            arcpy.Statistics_analysis(layer, path1, [["ID","COUNT"]], "ID")

            j = 0
            with arcpy.da.SearchCursor(path1, ["ID", "COUNT_ID"]) as cur:
                for row in cur:
                    if (len(row[0]) > 3) and int(row[1]) > 1:
                        print(row[0])
                        print(row[1])
                        arcpy.SelectLayerByAttribute_management(layer, "NEW_SELECTION", '"ID" = ' + "'" + row[0] + "'")
                        i = 1
                        with arcpy.da.UpdateCursor(layer.name, "ID") as outs:
                            for ids in outs:
                                while i < int(row[1]):
                                    ids = ("",)
                                    outs.updateRow(ids)
                                    print(i)
                                    j += 1
                                    i += 1
                            print("ID: '" + str(row[0]) + "' present in " + str(row[1]) + " and removed from " + str(i -1) + " features.")
            print("Total " + str(j) + "features corrected in " + str(layer.name) + ".")

But when I run the code, it just update only one feature with blank value and leave rest of them as it is.

What is a workaround for this?


texts: Can't install extensions nor grass-dev package

texts:
I've been trying to install several extensions on my newly installed GRASS GIS 7.2 (installed following instructions here), amongst them "i.fusion.hpf". The extension doesn't seem to matter as the problem is the same for everyone of them. So when I run g.extension, and try to add an extension, I get the message

ERROR: Please install GRASS development package

On g.extension's manual I read that this error meant I had to install the "grass-dev" package. I tried to do so but when I do it in the terminal, I get dependencies error messages... (google-translated from french)

~$ sudo apt install grass-dev
    Reading package lists ... Done
    Building the dependency tree
    Reading status information ...
    Some packages can not be installed. This can mean
    that you asked for the impossible, or, if you use
    the unstable distribution, that some packages have not yet
    been created or did not come out of entering.
    The following information will help you resolve the situation:

    Following packages contain unsatisfied dependencies:
    grass-dev: depends: libgdal-dev but won't be installed
               depends: libpq-dev but won't be installed
    E: impossible to fix problems, defective packages are in "keep as is" mode.

Here's what happens when I try to install libgdal-dev and libpq-dev:

$ sudo apt installs libgdal-dev libpq-dev
Reading package lists ... Done
Building the dependency tree
Reading status information ...
Some packages can not be installed. This can mean
that you asked for the impossible, or, if you use
the unstable distribution, that some packages have not yet
been created or did not come out of entering.
The following information will help you resolve the situation:

The following packages of unmet dependencies:
libpq-dev: Depends: libpq5 (= 9.5.8-0ubuntu0.16.04.1) but 9.5.9-0ubuntu0.16.04 must be installed
E: Unable to fix problems, packages are in "keep as is" mode.

I tried sudo apt install libpq5 but it's already installed.... strange. And sudo apt install -f didn't work out as well....

I tried uninstalling and reinstalling GRASS GIS but nothing changed... except now QGIS is uninstalled (so I reinstalled it, seems to work fine). I tried changing ppa. I went from ppa:ubuntugis/ubuntugis-unstable to ppa:grass/grass-stable and then downgraded to GRASS 7.0. Problem was even worse as I couldn't even access the add-ons database! Now I'm back on 7.2.

Starting to think it's because the packages are in development, somehow? Could the developers be updating them and thus they are inaccessible? I don't really know how that works to be honest...
Also on g.extension's manial/troubleshooting, it says:

  Since extensions have to be compiled on Unix based systems (Linux, Mac
  OSX etc.) unless a Python extension is installed, a full compiler
  environment must be present on the user's computer.

Does that mean I can install some Python extension? I'm really not sure what this means.


texts: FME error - failed to validate FME worker licence key

texts:
The FME licence I have access to is a floating one.

The issue seems to be easy to recognize (see attached log) however there are some doubts I would like to ask for.

While the process is really machine consuming I have some suspicious that the error of the FME is not due to "Failed to validate FME worker licence key" but it's more related to the task the machine has to deal with. Can it be ?

The error always occurs in the same point of proccessing time.

How does FME check for the licence key ? Just at the start of the workbench or there are some "milestones" during the processes where the licence has to be confirmed again ?

Please find below a part of the log.

Thanks in advance for the help.

Gregory

Failed to validate FME worker license key.
Program Terminating
Translation FAILED.
Failed to validate FME worker license key.
Program Terminating
Translation FAILED.
Failed to validate FME worker license key.
Program Terminating
Translation FAILED.
... Last line repeated 2 times ...
FME Session Duration: 4 hours 56 minutes 38.5 seconds. (CPU: 220.4s user, 11.6s system)
END - ProcessID: 1324, peak process memory usage: 1074456 kB, current process memory usage: 355076 kB
Error in child process '33', exited with code '2'
Program Terminating
Translation FAILED.
Warning: not all FMESessions that were created were destroyed before shutdown.  This may cause instability


texts: Remove the vertices of polylines which are too close to eachother

texts:
I want to remove the vertices or nodes of a polyline which are too close to each other. In other words, if the polylines are too close to each other, I want to cut both of them in the place where they are close to eachother ( see image).

I tryed v.clean from GRASS (break) but apparently it does not work for such a case.

First attempt

I tryed the the first recommanded methode (see below). but it does not work ar I wish. Hier I attach a new figure to show the details of the chosen parameter for v.clean.


texts:
I need to convert a raster into a polygon, and I ran into "Error 000863:  Invalid GP Data Type"
After searching, I saw this solution online.
Basically I needed to convert the float raster to an integer raster.  I used the times tool to multiply my raster by 1,000,000.0 to preserve all the decimal places.  When I compare the resulting raster, it seems Arcgis has rounded the value for each pixel up by one (See the pixel values in the images below). I was hoping to preserve these exact values.  Is there another approach I can try?


texts: how do i get a GDAL(using gdalwarp) clipped layer projected in pseudo mercator ESPG:3857 to not be shifted?

texts:
how do i get a GDAL rasterized layer projected in pseudo mercator ESPG:3857 to not be shifted ?
i have a vector layer projected in ESPG:3857,after rasterizing using this command:

gdal_rasterize -a_nodata -999 -te xMin yMin xMax yMax -a theField -tr 10 10 -l LayerName -co TFW=YES -co COMPRESS=LZW  VectorLayerPath outPutFilePath

the rasterization is perfect, but when i clip using this command:

gdalwarp -dstnodata -999 -q -cutline theMaskVectorLayerPath -crop_to_cutline -of GTiff inputrasterPath outPutRasterPath

the raster is shifted then clipped and the result is wrong.
is this happening to anybody else using pseudo mercator?


texts:
I am trying to open a geodatabase file and getting error message ValuError : Null Layer.
I am using geopandas.read_file.
import geopandas as gpd
gdp.read_file("AL_ALFL_conf_final.gdb")

I get error Null Layer: u''


texts: Detecting broken hyperlinks using ArcPy?

texts:
I am trying to use ArcGIS 10.2.2 to detect broken hyperlinks from a .shp and export them to a text file.

Is this possible to do this with Python scripting?

I only want to see the broken hyperlinks in the text file.

I tried using the script from Using ArcPy to verify paths of Hyperlinks?, but I could not get it to work.


texts: Is there any information in Sky Plot that is not in DOPs chart in terms of predicting accuracy in the field? (GPS data collection planning)

texts:
Joe only looks at DOPs chart and weather report for mission planning. Moe looks at both DOPs chart and Sky Plot and weather report. Is Joe missing out on something relative to Moe in terms of predicting accuracy in the field?

http://www.trimble.com/GNSSPlanningOnline/#/Settings

Note that there is a function to add obstructions or change cut-off elevation value in Settings, which affect both charts equally.


texts: Assigning values to polygon layer based on nearest point in point layer using ArcGIS for Desktop?

texts:
I have a point file depicting groundstations for rainfall measurment. I also have a polygon file representing districts of a state. I want to assign rainfall values in a particular district based on the station located there. I understand this can be achieved through a spatial join. However, I also want to assign these values to near by districts based on a value taken from the nearest station in the surrounding districts.

Please note that

a single district can theoretically have multiple stations.
I want to assign the exact value, and not an interpolated value.

I am working on ArcGIS Desktop 9.2.


texts: Selecting independent objects using GeoJSON on OpenLayers

texts:
I'm working with OpenLayers and loading a couple of GEOJson files which works perfectly. Now I want to be able to select individual components of the GEOJson file. As for example if I have more than one Feature (a couple of polygons) when I activate the hover selection, nothing happens.

I added the hover selector to this example http://openlayers.org/dev/examples/geojson.html

I can actually drag the elements independently but I can't select them. I assume it is because of the projection I'm using. Any Ideas how to fix it?

function init(){

var options = {
    projection: new OpenLayers.Projection("EPSG:900913"),
    displayProjection: new OpenLayers.Projection("EPSG:4326"),
    units: "m",
    numZoomLevels: 18,
};

map = new OpenLayers.Map('map', options );
var mapnik = new OpenLayers.Layer.OSM({layers: 'basic'});
map.addLayer(mapnik);

map.addControl(new OpenLayers.Control.MousePosition());

var lon = 12.18;
var lat = 53.81;
var lonLat = new OpenLayers.LonLat(lon,lat).transform(map.displayProjection,  map.projection);

var mapcent = new OpenLayers.LonLat(lon, lat).transform(map.displayProjection, map.projection);
zoom = 5.5;
map.setCenter(mapcent, zoom);
document.getElementById('noneToggle').checked = true;

var proj = new OpenLayers.Projection("EPSG:4326");

var file = [];
file[0] = "geojson/nuts3.geojson";

var style = [];

i=0;
while(i<file.length){
    style[i] = new OpenLayers.StyleMap({
        pointRadius: 2,
        strokeColor: 'red',
        strokeWidth: 1,
        strokeOpacity: 1,
        fillColor: 'red',
        fillOpacity: 0.2
    });

    var layer_name = "layer_" + i;
    vectorLayer[i] = new OpenLayers.Layer.Vector(layer_name,  {
        styleMap: style[i],
        projection: map.displayProjection,
        preFeatureInsert: function(feature){
        feature.geometry.transform(proj, map.getProjectionObject());
        },
    });
    map.addLayer(vectorLayer[i]);

    var gjson = new OpenLayers.Format.GeoJSON();
    var featurecollection = ReadFile(file[i]);
    //Readfile is just a normal Request.POST function
    vectorLayer[i].addFeatures(gjson.read(featurecollection));
    i=i+1;
}

controls = {
    selecthover: new OpenLayers.Control.SelectFeature(vectorLayer[0],{
        hover: true,
        highlightOnly: true
    })
};

for(var key in controls) {
    map.addControl(controls[key]);
    controls[key].activate();

}

            var report = function(e) {
            OpenLayers.Console.log(e.type, e.feature.id);
        };

        var highlightCtrl = new OpenLayers.Control.SelectFeature(vectors, {
            hover: true,
            highlightOnly: true,
            renderIntent: "temporary",
            eventListeners: {
                beforefeaturehighlighted: report,
                featurehighlighted: report,
                featureunhighlighted: report
            }
        });

        map.addControl(highlightCtrl);
        map.addControl(selectCtrl);
}

The geojson file looks something like

{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "id": 0, "properties": {}, "geometry": { "type": "Polygon", "coordinates": [ [ [ 9.244170, 48.847009 ], [ 9.246171, 48.822260 ], [ 9.248220, 48.796910 ], [ 9.248433, 48.794277 ], [ 9.290007, 48.791633 ], [ 9.298732, 48.785304 ],  [ 9.145165, 48.859655 ], [ 9.153553, 48.861241 ], [ 9.207493, 48.852769 ], [ 9.235667, 48.848344 ], [ 9.244170, 48.847009 ] ] ] } }
  ,
  { "type": "Feature", "id": 1, "properties": {}, "geometry": { "type": "Polygon",   "coordinates": [ [ [ 9.064848, 48.755577 ], [ 9.038994, 48.741317 ], [ 9.104720, 48.709080 ], [ 9.104721, 48.709079 ], [ 9.146561, 48.663574 ], [ 9.156822, 48.652414 ], [ 9.162166, 48.646602 ], [ 9.174173, 48.633544 ],  [ 9.064848, 48.755577 ] ] ] } }
,

... etc


texts: create sectors with given azimuths of one point in QGIS and show tabs info. following each sector direction

texts:
if there any way to create some sectors from one point with given azimuths,and show tab information on each sector following sector direction, or when i click the sector, all the information related to this sector can be shown.

the raw data:

ID  lon lat azimuths  sector_id

409970  113.243306  22.662132 9 10

409970  113.243306  22.662132 120 140

409970  113.243306  22.662132 256 260

the wanted style like this:


texts: I have an assignment to find out changed areas of the two different layers ( let's say 2000 and 2010).the layers are vector. Is it possible to detect the changed areas by using Arcgis?


texts: Installing QGIS on Mac OS 10.7. "The GDAL framework is required.", and yet I have it

texts:
what did I miss? I'm trying to install QGIS on my MacBook (Mac OS 10.7). I'm following these instructions: https://sites.google.com/site/eospansite/qgis-installation/installing-qgis-1-6-on-macosx

It says to first install GDAL following these instructions: https://sites.google.com/site/eospansite/installinggdalonmacosx

I installed GDAL but when I come to install QGIS I click on the QuantumGIS.pkg and it opens the install dialogue box but with a warning "The GDAL framework is required."

I know that GDAL is there, I can see it in /Library/Frameworks/GDAL.framework/ so I wonder if exporting the path didn't work. Anyone know how to check that? Or anyone know of anything else I may have overlooked? Thanks.


texts: Set raster layer additional no data value with PyQGIS

texts:
I need a raster layer transparency configured like this:

Since I load the layer with PyQGIS, the parameters should be set programmatically. For global transparency and transparency band I found that it works via myrasterlayer.renderer().setOpacity() and myrasterlayer.renderer().setAlphaBand(0), but I did not find a method to set the additional no data value to 0: http://qgis.org/api/classQgsRasterRenderer.html
How to set this value with PyQGIS?


texts: QGIS viewshed output is null -- incorrect

texts:
I am trying to run a very intensive viewshed in ArcGIS but university power outages are forcing me to use our SCC, which runs on Linux, so I need to use Grass and Q. I am running a much smaller test viewshed as a model for writing the code, and it's frustratingly incorrect. I know what the viewshed output should be -- it runs fine in ArcGIS and I've photographed from each of the OPs.

What I am looking for: a binary, cumulative viewshed with no outer radius, with each OP at a 1.65 m offset. (I realize the outer radius increases the computation, but it's necessary for my research).

What I get: a viewshed output that is completely null

I've checked the imported vector points, the DEM, re-run in ArcGIS, and even applied a small, 200 m radius. Nothing changes the output. I also don't get an error.

Can someone look at these parameters and check to see what I am doing wrong? I am very experienced with ArcGIS, but haven't needed to use Q or Grass very much.

(Apologies for my screen resolution -- GIS programs haven't caught up yet)


texts: How to set and configure a value relation edittype in QGIS 2.14?

texts:
I am trying to set the edit type of a vector layer to use a value relation, by pyqgis. So my question is exactly the same as this, except now it is qgis 2.14. I found that setEditForm is deprecated and this method is recommended instead:

vlayer.editFormConfig().setWidgetType(index,"ValueRelation")

setWidgetType sets the edit widget correctly but then I have the same problems as in the question regarding qgis 2.0. I.e. python commands are processed without errors or warnings but the when I check the properties in QGIS, the widget settings are empty:

vlayer.valueRelation(index).mLayer = target_layer.id()
vlayer.valueRelation(index).mKey = key_field
vlayer.valueRelation(index).mValue = value_field
vlayer.valueRelation(index).mAllowNull = True
vlayer.valueRelation(index).mOrderByValue = True

How do I define the value relation by pyqgis in qgis 2.14?


texts: Using contextual information to plan route?

texts:
I downloaded the .osm map data, and I imported it to postgresql using the m2pgrouting.
I have a CSV file on violence data, with three attributes, nature of occurrence, latitude and longitude.

How do I plan a route avoiding these points of violence?

What should I do, can I cross the data with street network?

Does pgrouting topology allow this?

If not, which approach should I follow?


texts: Finding local maxima with variable window search using ArcGIS Desktop?

texts:
I would like to find local maxima in different zones of raster file.

I got five zones (zone 1, 2,3,4,5) and would like to find local maxima with different window search for each zone in arcgis desktop 10.

For example, find local maxima with 3x3 moving window size  for zone 1, with 5x5 window size for zone 2, with 7x7 window size for zone 3.


texts: Can unicode characters be generated from unicode literals in QGIS expressions?

texts: I know that unicode characters can be copied and pasted into expressions in QGIS, but can they be generated from a unicode literal ie u'\u00C9' without using a custom Python function?


texts: Enabling a parameter in a field map in a Python script

texts:
I am creating a customized version of the Find Identical tool in Python. The script works just fine, but there's one thing that I can't seem to figure out.

In the standard ESRI Find Identical tool, the optional XY tolerance and Z tolerance parameters are disabled unless you choose the "Shape" field as the field you want to find identical values in. Then they become enabled. I'm trying to figure out how to do this in Python. I know this involves customizing behavior in the validation tab under properties, but all of the examples I've seen online show this for values and not something in a field map:

if self.params[3].value == "Get Spatial Weights From File":
    self.params[8].enabled = True
else:
    self.params[8].enabled = False

This code doesn't work for the field map. Is there a way to enable these parameters if "Shape" is selected in the field map?


texts: Convert point shapefile to polygon with QGIS

texts:
I am new to using QGIS. I have a point shapefile that I would like to convert to a polygon. I followed the directions for using point2one plugin but got this: .

This is what I want:

My attribute table has fields for ID, latitude,longitude, tracts, order(a field I created to group the IDs by tracts.


texts: Change Symbol Size and Color for Multiple Feature Classes

texts:
I need to expediently make an identical change to the symbol size and color of multiple feature classes. This can be done for layer transparency through a script such as:

mxd = arcpy.mapping.MapDocument("CURRENT")
for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.name == 'FC1' or 'FC2':
        if lyr.supports('TRANSPARENCY'):
            transparency = str(lyr.transparency)
            lyr.transparency = 45

But as far as I know, the color, size, and symbol type characteristics are not as easy to access and change. Any ideas on how this could be done?


texts: Can I change the precedence when combining "Control feature rendering order" with QGIS symbol levels?

texts:
I'm trying to produce a QGIS map style / symbology for processed openstreetmap data. For roads, I'd like to mimic how the OpenStreetMap web tile map "Standard" rendering style (produced with mapnik AFAIK) treats road intersections, road bridges over roads, road intersections on bridges and road bridges over road bridges etc.

Map tiles © OpenStreetMap contributors (CC BY-SA 2.0)

To produce the dark rim, I use the same trick as the "Standard" style: Draw a wide dark line and above it a bright line with slightly smaller width for the "interior" of the road. This looks alright for individual road features. To get continuation of roads (which might be split over several features in the data) and road intersections to look right, I can put all "interior" parts of these multi-level line symbols on a higher "Symbol level" (later rendering pass) than all the "rim" parts.

This is also suggested in the QGIS Training Manual. See their solution:

screenshot CC BY-SA 3.0 QGIS project

To get bridges on various levels drawn correctly over each other and over roads, I check "Control feature rendering order" and let it order by the z_order attribute as produced by osm2pgsql (cf. Applying OSM z-order after osm2pgsql with qgis). This works fine without the rim, but without the rim you can only really visually see the bridges when they are between roads of different type (which I symbolized in QGIS with different colors). For roads of the same type, the rim is needed to distinguish bridges from intersections.

Combining "Control feature rendering order" with "Symbol levels" gives a wrong-looking result, though. It seems that "Symbol levels" rendering passes take precedence over the feature rendering order: Bridges over roads start to look like intersections when all "interiors" are on a later rendering pass than all "rims":

OSM standard style for comparison:

Map tiles © OpenStreetMap contributors (CC BY-SA 2.0)

This can be fixed by putting both, bridge interiors and bridge rims on a later pass than non-bridge road interiors. But bridges over bridges will then still look wrong, as will intersections on bridges:

And here is the second example from above, rendered in QGIS:

If I use symbol levels with render pass numbers
non-bridge road edge < non-bridge road interior < bridge edge < bridge interior, it'll look like this:

I think they would look right if the feature rendering order took precedence over the "Symbol levels", i.e. order two strokes by rendering pass number only if their z_order value is the same, otherwise draw them in the order of the z_order. Can I get QGIS to do that? If so, how?


texts: Split vector polygon with a line - GRASS GIS

texts:
I was wondering if in GRASS GIS exists a way to split a polygonal area based on its intersection with a vector line, eg:

 ______                      ______
|      |                    |  *   |
|  *   |  + __________  ->  |______|
|      |                    |  *   |
|______|                    |______|

1 polygon  +  1 line    ->  2 polygons


texts:
Background: I am running QGIS 3.10.1, PostgreSQL 12.1 and PostGIS 3.0.0 on Windows 10. Using this environment, I am trying to run a simple query in DB Manager:

SELECT polygon('((10,20),(25,30),(30, 30),(10,20))')::geometry
UNION ALL
SELECT path('(1,21), (5,15), (9,20), (12,5)')::geometry
UNION ALL
SELECT box('(10, 21)'::point, '(16,10)'::point)::polygon::geometry
UNION ALL
SELECT circle('(20,10)'::point, 3)::polygon::geometry;

When I execute this query in DB Manager, it gives me the expected results:

I want to visualize the result of this query, so I try to save it as a layer, like so:

However, when I click the Load button, it gives me the following error:

When I look at the logs in QGIS:

Something is Missing

It is immediately apparent that the generated query is missing something inside the highlighted double quotes. What is the missing element? Is it the unique identifier, such as the _UID_ that I see in related posts on this same or similar topic?

QGIS Setting?

If so, is there a setting somewhere within QGIS where I can set this variable? Where is it?

I have seen multiple posts around this same issue, but none have pointed me to where I can set this variable correctly. It could be that version 3.10 is just too new for this.

An Experiment

As a little experiment, I created a table in the database with the following code:

CREATE TABLE ch02.geom_shapes AS (
  SELECT polygon('((10,20),(25,30),(30, 30),(10,20))')::geometry
  UNION ALL
  SELECT path('(1,21), (5,15), (9,20), (12,5)')::geometry
  UNION ALL
  SELECT box('(10, 21)'::point, '(16,10)'::point)::polygon::geometry
  UNION ALL
  SELECT circle('(20,10)'::point, 3)::polygon::geometry
);

When I open this table directly in QGIS, I get the following:

No errors!

This is precisely the result that I was expecting to see. I do hope that someone knows how to get around the error that occurs when trying to create the layer from DB Manager.


texts: Point.in.polygon function in sf package miscalculate the distribution of points

texts:
I am using point.in.polygon function in sf package to determine whether a set of data fall within a polygon or not.
Here is a snippet of the data and polygon visualization

As the figure shows, there must be at least one point fully lies in the polygon, however, the point.in.polygon function returns all zero, which means the function thinks all points are exterior to the polygon.

Where am I do wrong?

codes
subDist <- point.in.polygon(CD8_Pts[,1], CD8_Pts[,2], sub_current_poly[,1], sub_current_poly[,2])

Edit: The following link is the polygon coordinates.
https://livejohnshopkins-my.sharepoint.com/:u:/g/personal/hmi1_jh_edu/EeC-sUzqPZJFtoUzM-r5DQ4B_aIsLri4G29RAi47uO42zw?e=Gw2oHt


texts: Import JSON or CSV with corner coordinates

texts:
I am trying to import meteorological forecast data into QGIS. Unfortunately the data lacks coordinate attributes, the spatial data is specified by the met institute as corner coordinates (SE, SW, NE, NW).

The data is in .JSON-format, but I managed to convert it into .csv as well. I tried to import as comma delimited text file but failed to add the georeferencing.

Any ideas?

Here's a link to the data if it helps:

http://opendata-download-metfcst.smhi.se/api/category/pmp1.5g/version/2/geotype/point/lon/16/lat/58/data.json


texts: Legend properties for a raster: how to get rid off one value?

texts:
I have a raster in ArcMap, which has the values 0,1 (it was derived from Raster Calculator)
I want to put in the Legend only the value 1. I put no color and no line to value 0, but I have this big gap between the Layer name and the label (see the following picture).
How can I get rid off the value 0, in order to have a nice looking legend?


texts: Adding a Vertex To Exisiting Feature(linear) Using OGR Python

texts: I am working on a project that involves finding the intersection of linears in a shapefile, if that vertex(the found intersection) doesn't already exist for each(or one) of the features, then that vertex(the found intersection point) should be added to both features. I'm using the shapely and OGR modules for python. Can somebody just give me a general example of adding a vertex to an existing feature in shapefile using either or and saving it. Thanks.


texts: Pick point X,Y values getting from ArcMap window

texts:
How to pick a point on Arcmap window using with ArcObjects(BaseCommand)? I want to get X,Y values of a point. I have tried to pick point using with below code. but debugging time not showing arcmap window form picking point.

public ESRI.ArcGIS.Geometry.IPoint GetScreenCoordinatesFromMapCoorindates(ESRI.ArcGIS.Geometry.IPoint mapPoint, ESRI.ArcGIS.Carto.IActiveView activeView)
        {
            if (mapPoint == null || mapPoint.IsEmpty || activeView == null)
            {
                return null;
            }
        ESRI.ArcGIS.Display.IScreenDisplay screenDisplay = activeView.ScreenDisplay;
        ESRI.ArcGIS.Display.IDisplayTransformation displayTransformation = screenDisplay.DisplayTransformation;

        System.Int32 x;
        System.Int32 y;
        displayTransformation.FromMapPoint(mapPoint, out x, out y);
        ESRI.ArcGIS.Geometry.IPoint returnPoint = new ESRI.ArcGIS.Geometry.PointClass();
        returnPoint.PutCoords(x, y);
        return returnPoint;
    }


texts: Cropping Landsat images with shapefile in MATLAB

texts:
I want to crop several Landsat images, of different satellites, in MATLAB employing a shapefile.
For only one image, I could manage to do it properly.

But if I try the same with a second image, the matrices associated to the images are of different size, mainly in the black part (zeros) while the mask (white part) in general keeps the same size. Also the location of the mask in both of the images differ within the figure.

All in all, if images have matrices that do not have same number of pixels and the masks are not located in the same part of the figure, it make hard to compare and make operations with them.
So I was wondering if there is a way to improve the code above to enable it for cropping several images with different extents using the same shapefile, while keeping the mask in the same location for all of the images as well as the same total number of pixels including the zeros of the black part around the image.

I1 = geotiffinfo('MyImageL8.tif');

[x,y] = pixcenters(I);

[X,Y] = meshgrid(x,y);

roi=shaperead('MyShape.shp');

figure('color','w');

mapshow(roi);

rx=roi.X(1:end-1);

ry=roi.Y(1:end-1);

mask = inpolygon(X,Y,rx,ry);

figure;

imshow(mask);

G1 = imread('MyImageL8.tif');

G2 = imread('MyImageL5.tif');

figure;

imshow(G1);

imshow(G2);

CUT_G1 = im2double(G1).* mask;

CUT_G2 = im2double(G2).* mask;

figure;

subplot(1,2,1)

imshow(CUT_G1);

subplot(1,2,2)

imshow(CUT_G2);


texts: ArcGIS Esri feature layer query get object lat and long

texts:
I am trying to do a query on a feature Layer and the get an object's Lat and Lng, so I can center it on the map. However after doing the query, I get a graphic that has only an x,y position. the extent.center.lat/long is null

Here is the code:

    centerAndZoomOnAsset(Id: string) {
    let query = this.pipeFL.createQuery();
    query.where = `Id = '${Id}'`;
    query.outFields = ['*'];
    this.pipeFL.queryFeatures(query).then((result: esri.FeatureSet) => {
      debugger;
      const foundFeatureGraphic: esri.Graphic = result.features[0];
      if (foundFeatureGraphic) {
        const center = foundFeatureGraphic.geometry.extent.center.clone();
        //this.mapView.goTo(foundFeatureGraphic);
        this.mapView.goTo({target: center});
      }
    });
  }

center.latitude and center.longitude are both null.


texts:
The plugin installs fine but fails to load with the following Python error

cannot import name QgsImageOperation


texts: How to simplify osm2po's routable graph?

texts:
I'm building a pgrouting routable network thanks to osm2po, the problem is that my network is too big. Indeed I want to create isochron and 1 node = 1 canvas circle, so I don't want too many nodes.

I managed to reduce my network with grass (v.clean) but when I do that I'm losing all my costs (and reverse_costs).

Is it another solution?


texts: GeoServer PointSymbolizer graphic artifact when zoomed

texts:
I'm using GeoServer 2.14.0, I created a red circle using PointSymbolizer with size=35. But in some cases, the map shows this artifact when I zoom in:

Basically the map still shows parts of the circle rendered with the previous zoom. I guess the problem is the tile render/refresh, a problem similar to Leaflet.js rendering 'chopped' WMS point symbols for variable sized points or http://osgeo-org.1560.x6.nabble.com/Geoserver-cut-off-symbols-close-to-tile-edges-td3798283.html

Any changes to the GeoWebCache settings (metatiles size and gutter) didn’t solve the problem. After a few tries, I realized that the larger the point size, the more evident the artifact is.

Do you have any suggestions?