texts:
I have been trying to learn to work with QGIS (open-source GIS software) through some tutorials available online for classification of image for land cover. These tutorials used a dataset from Landsat 8. The tutorial used the semi-automatic plugin v 3.0

Now I want to use the same process for an image in .jpeg format that I downloaded from internet.
Step one involves converting the raster bands from DN to Reflectance. For this a .txt file (a metadata file named .MTL) is required. How do I generate this .txt file for my image.


texts:
I am trying to transform a local vector rECEF from ECEF coordinate system to ENU coordinate system (rENU) where

ECEF = earth-centered earth-fixed coordinate system

ENU = East, North, Up coordinates, a local Earth based coordinate system

I saw in more then one book the following equation:

where

,

and

What i don't understand is why

since  and

What am i missing?


texts:
I have a collection of polygons in a geojson format and I'm trying to find for one polygon, its neighbors.
I tried to use the $geoIntersects operator like so
coords = [...] # some coordinates of a polygon of my collection
db.coll.find({
  "loc": {
    "$geoIntersects": {
      "$geometry": {
        "type": "Polygon" ,
        "coordinates": coords
      }
    }
  }
})

It didn't give me any results.
For clarification here two polygons that are next to one another and I'm trying to find one when I'm calling my function with the other.
They do not need to share exactly one full side, just a portion of it.
Is that something that's doable with mongo? Or should I aim for another strategy?


texts: Create polygons from polylines with polylines attributes

texts:
I am working in ArcMap 10.3.1 and I have a feature class with some 500,000 polylines, from what I need to generate polygons, that have to keep attributes from the polylines. I have tried the "Construct Polygons" tool in Edit mode, but it returns empty attributes table, and "Feature to Polygon", but it doesn't return attributes at all. I also tried the spatial join for the generated polygons with the lines feature class, but there are too many errors with attributes being assigned to incorrect features.

Is there any other method I can  use to achieve the goal of keeping attributes in the generated polygons?

Edit: added a couple of snapshots of line feature class fragment and, yes, they do overlap.


texts: Extracting specific feature class from dataset in ModelBuilder?

texts: I would like to be able to run a model that iterates through datasets and at one point extracts just a single point feature from within the dataset and save it in a different geodatabase.  The problem I run into is that the feature class to feature class tool will not recognize my inline string for the path of the feature class.  It wants me to point a feature class within the model to the feature class to feature class tool.  Is there a way to "extract" a single feature class from a dataset and then use it as the input for the next tool in the chain within model builder?


texts: Convert pixel values of a raster from meters to feet and from feet to meters - GDAL

texts:
I want to convert the elevation of a GeoTIFF from meter to the feet, so I can create a contour in feet. I've used gdal_translate:

gdal_translate -scale 0 0.3048 0 1 inputmeters.tif outputfeet.tif

but when I create a contour, both are identical. and with the gdal_info I don't see any changes in the metadata.
I'm not sure if the gdal_translate is the command that I need!

It is necessary for me somehow to check if the outputfeet.tif unit type is feet.


texts: Making segmented roads appear as one line in QGIS?

texts:
I am footling around with QGIS and the British OS Strategi dataset. If I render the roads with a simple line style they look fine:

However, if I apply something slightly fancier, it looks wrong, because the style is applied to each segment of the road as opposed to the whole road (as we perceive it), as seen in the white lines on the blue motorway below:

How can I merge these segments together, or failing that force QGIS to treat them as one line?

I'm not using the data for navigation so I don't mind if the segments are lost.


texts: Converting a .lasd (Lidar) coordinate system?

texts:
How can one change a coordinate system for a .lasd (Lidar) dataset?

The file is currently in UTM_Zone_15N and I am trying to convert it to IL_East.


texts: PostgreSQL - Merge only lines that share a specific attribute

texts:
I have a trail that is split up into multiple small parts. On some of the parts i have managed to get the traile name as attribute "text". The rest has text as "null". Now I would like to join all lines in the trail that has close enough endpoints and that has either share the "text" or has the text as null.
Example of attribute "text" in lines that all have close endpoints:

null
null
Trail1
null
Trail1
null
null
Trail2
null
Trail2

In the above list I would like to group:

1-7 with a common attribute "text" = Trail1
8-10 with a common attribute "text" = Trail2

EDIT - Added illustration(Left part of the long black line should be grouped as Trail1, right part as Trail2:

I have used ST_ClusterDBSCAN(geom, eps := 5, minpoints := 2) to group by distance before but not on attributes like this. If i would run it on the example above it would put all lines in the same entity Wich is an issue.
The queryes i have been using this far.
1.This query joins two tables, one with all of the small trail parts and one with the texts, so im joining the texts that are within 50 meters from a trail and assume that is the trail name. This gives me the data in the illustration above.
DROP TABLE IF EXISTS "lines_with_text";
CREATE TABLE "lines_with_text" AS
SELECT
  lines.id,
  max(lines.location) as location,
  max(texts.text) as text,
  lines.geom
FROM
  "line_table" as lines
LEFT OUTER JOIN (SELECT * FROM "text_table") AS texts ON ST_DWithin(texts.geom, lines.geom,50)
GROUP BY lines.id, lines.geom

2.Then i tried to merge them with ST_ClusterDBSCAN but this does not give me the desiered result. This will join all lines that are close together an will split them when two lines are not within 20 meters. I would instead split when a new trailname is found rather than when the distance is to long. Usally the lines are between 0-2 meters apart from each other(they are not perfectly aligned but sometimes they are more than 20 meters apart and then i have accepted that it will be a new polyline. Maby with the same Trail name as the prior one.
DROP TABLE IF EXISTS with_text_joined;
CREATE TABLE with_text_joined AS
SELECT
  min(location) as location,
  ST_Union(geom) as geom,
  min(text) as text
FROM (
  SELECT
    location,
    geom,
    text,
    ST_ClusterDBSCAN(geom, eps := 20, minpoints := 2) OVER() AS cluster_id
  FROM "lines_with_text" as lines
) sq
GROUP BY cluster_id


texts: Changing Text Element Leading Causes Position Changes

texts:
I have several maps I need to update that have a long list of text contained in a single text element. I wanted to vertically compress the text, while retaining the font size. The problem is if the leading is set to anything other than 0, the entire text body moves position on its own whenever I pan around in layout view.

Is there a setting I am missing that could be causing this? Or is this a known bug?


texts: How can I upload 10 million points with 100,000 categories?

texts:
How can I upload 10 million points with 100,000 categories and have CartoDB randomly assign a color to each cateogry?

Currently, it only seems to color the first 10. And from what I've seen on other questions the only way to color more is to edit the CSS files. I don't think this is too feasible for what I'm trying to do. Is there anything built in for this use case?

The reason I have so many categories is because these points cover all of north america and are clustered together based on their category. So for example the items can be households (one points per household) and the category would be the neighborhood the household is in. I actually don't care about having unique colors, so long as two nearby clusters don't share the same color.

I decided to just add a new random category value to determine the coloring and add my own CSS for the ~50 unique colors I chose.


texts: Delimit Hydrographic Basin, based on a single river

texts:
I find myself doing a hydrological study, where I have to delimit the river basin based on a single river (polyline).

I have the DEM where the river of my interest is, and I have the river as a shapefile, but generally when I have made the delimitation of basins in ArcGIS, I have done it based on a network of drainage, but not on the basis of a only river of interest.

Is there a tool in which you can do all the hydrological processing, to obtain as a result the river basin of a specific river?


texts: Sharing ArcMap with router and simultaneously two or more than two users can data entry in attribute table?

texts: I have two computer, one router and I want to share ArcMap and simultaneously from two computer filling the attribute table.


texts:
I don't have Pathfinder on my computer. I need to pull data off of my Trimble unit (Geoexplorer 6000 series) and save it on an external hard drive or USB key, without using Pathfinder. This way, I can send the storage device to a colleague in another location who has Pathdfinder and can import it properly and convert it to a shapefile for me. Is this possible? I've tried to copy and past raw data from the unit to my computer and have been met with an error.

Sending the GPS unit back and forth to import data from it is quite expensive and I'm trying to avoid this.


texts: Copying geometry from one feature to another feature with NULL geometry using QGIS?

texts:
I have a shapefile in which some of the features have NULL geometries. What I would like to do is to copy existing geometries from other features and paste them to those features within the same layer that have a NULL geometry.

Is there an easy way to do this within QGIS, preferably the GUI?


texts: Raster manipulation in ArcGIS - Flood Mapping with HAND

texts:
I am trying to create flood maps using the height above nearest drainage (HAND) method using data from the continental flood inundation mapping (CFIM) project. I have the following data:
A raster map of subcatchments

A map of HAND elevations

A table of flow/water height for each catchment

I will calculate a flow for each catchment, then using the table, I need to apply a unique water height to
each subcatchment. I am new to ArcGIS and I am looking for general guidance for how to do this. I need to break up the map of HAND elevations into the subcatchments shown, then apply a unique elevation mask to each layer.
How should I do this?
I can code in Python but I'm not sure what functions I need to use.


texts: Balancing multiple supply and demand locations in QGIS

texts:
I am looking for a way to balance supply and demand with respect to distance, between two vector layers in QGIS.

I have two vector layers 50x50m grid

Supply (The red negative values)
Demand (The green positive values)

The plan is to develop a map to show arrows with values to indicate where all the positive (green) demand values need to acquire the negative (red) values from supply to =0.

Some Additional Notes:

Develop the shortest routes possible
Cells can be split multiple directions
Cells can receive units from multiple directions

Question: If I want to perform this Network Analysis in QGIS what is the proper terminology for it?


texts: Can't save result of "Extract Values to Points (Spatial Analyst)" inside the working feature dataset!

texts:
I have a file geodatabase in which I have a feature dataset. I have some feature classes inside the feature dataset. I wonder why when I attempt to save the result of Extract Values to Points (Spatial Analyst) inside my working feature dataset, the result is genereated inside the file Geodatabase but not under my feature dataset. No matter what I do I can't make it save the result inside my feature dataset. I can easily save the results of for instance Extract Values to Points (Spatial Analyst)  inside my working feature dataset. But "Extract Values to Points (Spatial Analyst)" seems to be different in that regard. Any idea why this is happening?

Edits:
I am using Arcmap 10.0 calling this tool from Arc toolbox. All the layers are in the same coordinate system. Here is a snapshot before I run the tool:

But as can be seen here the result is generated outside the feature dataset:


texts: Which label from different points inside a polygon is the most repeated

texts: I have millions of points, each point with a label. These points are inside of different buffers (polygons). I want to create a new column that puts the value of the label that is repeated more time inside of a polygon. I draw something just to give you an idea. The circles are the buffers, and the labels of each point are 1,2,3 in the buffer on top left I gave it the label 2 because we see two points with that label.


texts: How to spatially join several attributes of one layer to a unique one on a second layer

texts: My goal is to link all the attributes of an input layer (that has a UTM CRS projection matrix and point coordinates but no geometries) to the geometries of another layer based on the spatial analysis. In particular, I want to join them if they intersect or if the geometry contains the dot equal to the coordinates in the input layer. Using the plugin MMQGIS, I am able to merge the two dataset. The problem is following: often multiple observations lies on the same geometry. MMQGIS only returns a random attribute of the input layer but reports correctly the number of points that intersect the geometry of the second one. I would like to have instead a new layer in which all the input attributes are reported. This means that the geometries' id will be duplicated in those cases with multiple inputs on the same geometry. Is there a way to easily perform this task? Both PyQGIS or icons solutions are useful.


texts: Heat map displaying as one solid layer with no difference based on point density

texts:
I'm trying to create a heat map in ArcMap 10.4.1. I've already selected a base map, created a shapefile out of my data points, and overlaid that shape onto the map. Next I go tools -> point density. I set Input Point Features to my shapefile. I click on environments, then Raster Analysis, and then I set mask to my shapefile. Finally I go to Processing Extent and set the boundaries I want using another layer.

After all of this, I get the entire area filled in with one solid uniform color. The densities of the points are not reflected at all. Also, when I look at my table of contents, I see that the ranges for every single color are the same.

Can anyone tell me what I might be doing wrong?


texts: How to get map from GeoServer using longitude, latitude, zoom level, device width and height?

texts:
I need show a map as image (gif) in my Windows mobile. Only I will get center lat/lon value, zoom level and  width-height mobile screen.
Please help me to solve this.


texts: Error when Projecting Rasters using ArcPy?

texts:
Using ArcGIS 10.3 I'm receiving the dreaded 'ERROR 999999' when attempting to project a geodatabase full of rasters from WGS84_UTM_zone_10N to Nad83_BC_Environment_Albers

import arcpy

InRasterGDB = r"C:\Users\PoirierM\Documents\GISData\LandsatSSTImagery.gdb"
OutRasterGDB = r"C:\Users\PoirierM\Documents\GISData\LandsatSSTImageryAlbers.gdb"
OutCoordSystem = "3005"
ResampleType = "BILINEAR"
CellSize = "30 30"
GeoTransform = "WGS_1984_(ITRF00)_To_NAD_1983"
InCoordSystem = "32610"

arcpy.env.workspace = InRasterGDB

for fname in arcpy.ListRasters():
    print fname
    print OutRasterGDB + "\\"  + fname
    arcpy.ProjectRaster_management(InRasterGDB + "\\" + fname, OutRasterGDB + "\\" + fname, OutCoordSystem, ResampleType, CellSize, GeoTransform, "", InCoordSystem)

This should work, but it returns the following error when run in the Python Window in ArcCatalog:

Runtime error
  Traceback (most recent call last):
   File "<string>", line 16, in <module>
   File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", line 9009, in ProjectRaster
     raise e ExecuteError: ERROR 999999: Error executing function. Failed to execute (ProjectRaster).

I can do a single raster manually through either the tool interface or within modelbuilder. Can anyone suggest why it won't work as a Python script?

I've got a working script now - see my (MPoirier) answer below.


texts: Exporting multiple shapefiles (mixed geometry types) into one KML?

texts:
Is there a way, using open source software, to export multiple shapefiles of mixed geometry types, and different attributes, into one KMZ/KML file?

Something similar to from ArcMap "Map to KML", but styling information is not necessary.


texts: How do I create a target geometry to use with JTS GeometrySnapper?

texts:
I am dealing with geospatial queries on a database. The query returns all data sets within an arbitrary bounding box, given a point associated with the data set is within this bounding box.

Since the bounding boxes provided have very high precision, I was not able to apply caching mechanisms based on the bounding boxes latitude and longitude values. In order to be able to cache queries, I thought it would be a smart idea to normalize bounding boxes before querying the database.

My initial idea was to simply truncate lat/lon decimal places, but this approach doesn't allow for configuration of the precision.

Later, I found JTS GeometrySnapper which allows to snap to arbitrary geometries.

I would love to create a geometry that allows me to snap to a grid comprising of lines every n meters both in latitude and longitude. How do I construct such a geometry?

This is my approach so far:

  private def sanitizeBoundingBox(boundaries: BoundingBox): Geometry = {
    val precisionModel: PrecisionModel = new PrecisionModel()
    val geometryFactory: GeometryFactory = new GeometryFactory(precisionModel, 4326)
    val box: Geometry = geometryFactory.toGeometry(new Envelope(boundaries.minLongitude, boundaries.maxLongitude, boundaries.minLatitude, boundaries.maxLatitude))
    val snapper = new GeometrySnapper(box)
    snapper.snapTo(???, ???)
  }


texts: How can I collect geometries from a singlepart shapefile to a multipart using GeoPandas or other non-ArcPy solutions?

texts:
I've tried to use the geopandas.collect.tool to do this using the snippet below, but it clears the fields from the input shapefiles and does not successfully collect geometries.
doc: https://geopandas.org/reference/geopandas.tools.collect.html
import geopandas as gp

gdf = gp.read_file('output/input.shp')
geos = gp.GeoSeries(gdf)
gp.tools.collect(geos)
gp.GeoSeries.to_file(geos, 'output.shp')

Can anyone correct my method or suggest a better solution?


texts:
I'm sorry for my English.
I look for tool or script for combining text fields in this way:
e.g.

field_1 (contains 3 to 6 signs)

field_2 (contains 6 to 10 signs)

[field_1] & x & [field_2]

where x = variable number of signs (e.g. "-" or "0") which will give me globally always 25 characters as a result of combining.
I'm working in arcgis or qgis. It could be also tool in acces or excel.


texts: Calculating coordinates of points from just distances, no bearings

texts:
The task I'm giving myself is about enabling people with little equipment to produce a map of a forest, or a garden with an important presence of trees.

I participate to openstreetmap and I plan to (enable people to) contribute tree locations for botanic gardens in areas where Bing aerial coverage is not available, or where it is useless because tree canopy is too thick to allow for identification of individual trees.

My hope is that there is a QGIS plug-in (and if it does not exist, to develop one myself) that accepts an amount of reference points distributed around the area of interest, a matrix of distances among points (referenced or not), and populates a vector(points) layer with the calculated 3D coordinates.

I have been looking at trilateration, and I found this 22 years old article.

I suppose the math still holds.

I do not want to work with angles, only distances.

not yet an answer, but I have been working at it, and I have a couple more details about the issue.

I need my input coordinates in metres (or feet, or millimetres), not in degrees, because I need an isometric coordinate system.
The math in that article still holds, and numpy helps a lot, specifically, the article suggests solving the overdefined problem by least squares, and numpy contains the function np.linalg.lstsq.
I gave up hope to work in 3D and I am going to assume earth is flat, but really flat like the Netherlands, or at least that the digital device I will acquire to measure distances is able to split distances in their horizontal and vertical components.
I have a working prototype, but I still need to discover how to create a point feature using the coordinates in the coordinate system of the layer.

the code for the coordinates calculation is not at all long, while the difficult part was about choosing which point I can reference given the points I have already calculated:

def find_point_coordinates(points, distances, point_id):
    import numpy as np
    from numpy.core.umath_tests import matrix_multiply

    connected_to = [id for id in sorted(distances[point_id]) if points[id].get('coordinates')]
    connected_matrix = np.array([points[id]['coordinates'] for id in connected_to])
    A = connected_matrix[1:,] - connected_matrix[0,]
    A = 1.0 * A  # make sure we work with floating point values

    ## squared distances vector, beacon_i to first beacon for which we have distances
    D_i1_2 = matrix_multiply(A*A, [[1],[1]])
    ## distances of targeted point from used reference points
    dfb_sel = np.array([distances[point_id][ref_id] for ref_id in connected_to])
    r2 = dfb_sel * dfb_sel

    rhs = ((r2[0] - r2[1:]).reshape(D_i1_2.shape) + D_i1_2) / 2.0
    r1, r2, r3, r4 = np.linalg.lstsq(A, rhs.reshape(rhs.shape[:1]))
    return connected_matrix[0,] + r1


texts: Duplicating active layer and adding definition query using PyQGIS?

texts:
Using: QGIS version 2.12.0-Lyon. No Add-on available.

Problem: I do not know where to start after opening the Python console with regards to writing a script that creates a new layer per value queried from a field within QGIS Desktop. I have experience with arcpy and am trying to start with a basic scripted action to learn Python for QGIS.

At the moment I can list all the attributes found in the "active" layer see below:

    lyr = iface.activeLayer()

    features = lyr.getFeatures()

    for feat in features:
        attrs = feat.attributes()
        print attrs[1]

However, I would like to iterate this and create a new memory layer (i.e. only available whilst the project is open) per feature per attribute creating a new layer where a query is set to equal that attribute (i.e. Name = feat[0])

Can anyone assist in building a template script solution?

This can be partly completed without using Python by using the "Split Vector Layer" tool. However, this tool will create a folder of shapefiles and not "in memory" layers.


texts: Grouping polygons that are part of fishnet using ArcGIS Desktop?

texts: If I have the grid below, which consists of 4,000 squares, how could I aggregate the squares such that I end up having 1000 larger squares, each of which contains 4 of the original smaller squares, without overlapping?


texts: Calculating difference between two rasters using ArcGIS for Desktop?

texts: I'm using ESRI Desktop 10.2 with a spatial analyst extension. Can Spatial Analyst be used to calculate the difference between two water level raster files, or do you have to use 3D Analyst?


texts:
I'm trying to create dynamic leader lines by using a PostGIS view in addition to the QGIS „Move Label“ tool.

CREATE VIEW leader_line AS
SELECT
gid,
ST_MakeLine(geom, ST_SetSRID(ST_MakePoint(xcord_label, ycord_label), SRID))::geometry(linestring, SRID) AS geom
FROM point
WHERE xcord_label IS NOT NULL;

This works fine for all labels WHERE ST_X(geom) < xcord_label but creates wrong looking leader lines for labels WHERE ST_X(geom) > xcord_label.

Does anybody know how to get properly placed leader lines for labels WHERE ST_X(geom) > xcord_label? Is there any way to refer to the xmax coordinate of the labels?


texts: how to delete a style that associated to a layer?

texts: I want to apply a style to a layer by rest.at the first time everything is ok. But when I change the content of the sld file in the style folder mannually,it doesn't apply to the layer. I want a way to do  this by rest,not by submit the style through the style editor in geoserver.and How can I remove binding of astyle from a layer?


texts: Making meter and Arcmin based raster compatible for analysis

texts:
I am trying to work with temperature and precipitation data from worldclim (Bioclim) and the Corine Land Cover data. however, the rasters are in two different grid sizes and formats, 2.5 Arc Min and 100m respectively. I need to know how to convert one grid format into the other or some way of analysing using both rasters if possible.I know arc mins are not exactly a straight forward size and it is dependent on the latitude so I am not sure what to do.

I am using ArcGIS 10.3 for Desktop.


texts:
How do I use python to find raster layer properties like
The GDAL driver and/or description used to load a raster in QGIS
For Example : for a Tif 'GTiff' and 'GeoTiff' or
       for an ESRI Grid 'AIG' and 'Arc/Info Binary Grid'
For a vectors you use layer.dataProvider().storageType() but when applied to a raster layer it returns AttributeError: 'QgsRasterDataProvider' object has no attribute 'storageType'
I'm using QGIS 3.10


texts: Using SCP image classification to find swimming pools?

texts:
i'm currently trying to use the scp image classification plugin to find swimming pools in an aerial image. It is just an RGB image as if it had infra red this'd be a lot easier.

Does anyone know how i could improve my results in terms of image classification without having access to imagery with IR or NIR bands using this or another classification plug-in for QGIS?

The three algorithms used by the scp plugin were unable to provide raster files demarcating the swimming pools using RGB imagery.


texts: How do I limit update/insert/delete capabilities for a WFS layer in a QGIS project?

texts:
In QGIS 2.0 I have a layer based on a query to a PostgreSQL/PostGIS database and would like to serve this as a WFS layer. Users of the layer should be able to select subsets of the data based on properties such as dates or unique IDs, but should not be able to update, insert, or delete records.

The Creating a WMS/WFS from a QGIS project section of the documentation simply says, (emphasis mine)

  In the WFS capabilities area you can select the layers that you want
  to provide as WFS, and specify if they will allow the update, insert
  and delete operations.

How, exactly, do I specify this? What is the correct way to deny, and to allow, these actions?

Here's what I've done so far. I added the layer I want to publish to a new QGIS project. Then in the Project Properties window, under OWS server, I checked the box in the Published column for the layer. Now I don't know what to put in the free-form text input boxes in the Update, Insert, and Delete columns. I've tried putting values of "0", "NO", and "FALSE" in here, as well as more facetious entries such as "Sure" and "Why not?". No matter what I enter though, when I click Apply or OK, QGIS gives no error or feedback of any sort. When I return to the project properties later, the text boxes are all blank as they were initially.

Oh, and in case it makes a difference, I'm running this under Windows (currently 7, but it will be deployed on Server 2008 or 2012 eventually).


texts:
Unfortunately, GDAL does not support a particular projection I am trying to achieve (the ob_tran rotated pole projection), so I am trying to use PROJ4 to convert a geotiff file to the projection I want.

However, I am not sure what the correct syntax for such an operation is. Using the command:

invproj -m 57.2957795130823 +proj=ob_tran +o_proj=latlon +o_lon_p=20 +o_lat_p=20 +lon_0=180 input.tif

I get only a lot of gibberish in the console, such as:

133d13'9.044"E  62d0'32.793"N  <Item name="missing_value" sample="0">-32767</Item>
133d13'9.044"E  62d0'32.793"N  <Item name="NETCDF_VARNAME" sample="0">Band1</Item>
133d13'9.044"E  62d0'32.793"N</GDALMetadata>
133d13'9.044"E  62d0'32.793"N
133d13'9.044"E  62d0'32.793"N� �D�k�����������k�R�f�����P�w��B�����|�[�=�$

This is not due to the particular projection I am trying to use. Trying a normal Mercator projection with the command:

proj +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +units=m +k=1.0 +nadgrids=@null +no_defs input.tif

I also get a lot of gibberish.

0.00    -0.00�4�L�}�z���`����)�,�$��������
0.00    -0.00��6�������G�R�f�@�
0.00    -0.00
0.00    -0.00�
0.00    -0.00���������3���<�����;
0.00    -0.00�����

What am I doing wrong? From looking at the manual, I thought I need to put the name of the file(s) at the end of the command line, but apparently, that is wrong?

How do I convert a file to another projection with PROJ4?


texts: Comparing two point clouds, one in WGS84 and the other in ECEF

texts:
I'm really a newbie and this task seems over my head. Suppose you have one cloud from LiDAR (a LAS file) with coordinates in WGS84 and another cloud, obtained from structure-from-motion, in ECEF. They represent a 100 mt long street and one (say the first) have to be used as ground truth for a statistical comparison.
My intention is to use the Python library laspy to extract the coordinates of the points and pyproj (the Python interface to PROJ) to convert from WGS84 to ECEF. I'm especially worried about losing accuracy when doing this conversion. One example of the coordinates I have in my LAS file (1.2 version, point data record 1, proj=latlon) is:
((-8495308, 1344287, 84460, 100, 9, 1, 0, 0, 1576, 1.52689438e+08),)

with scale [1e-08, 1e-08, 0.0001]
and offset [30.40828783662081, 59.9284293900385, 20.370796515028587].
In order to convert to ECEF, I have to apply both scale and offset to obtain actual values of lon, lat and elevation. Am I losing accuracy doing these operations? Is there any alternative?


texts: In ModelBuilder use Parse Path (path and name) to get output FC name

texts:
I built a model (in ArcGIS 10.6) and the Parse Path tool is working well if I using the Name type.

But I need to use the input .gdb name also into the output feature class name. e.g.

Inputs:

C:\clip\seeq.gdb\v1; C:\clip\lake.gdb\v1; C:\clip\lake.gdb\v2; C:\clip\lake.gdb\v3; C:\clip4q\hill.gdb\v1...

and I like to get these output FC names:

seeq_v1; lake_v1; lake_v2; lake_v3; hill_v1...

How could I use path? Or could you know something helpful?


texts: How can I clean/fix polygons in a Node.js application?

texts:
I am working on a simple Node.js application where users can walk around an area with a mobile device using GPS. This track is stored as a Polygon, but sometimes (due to inaccurate GPS signals or users not walking strictly around something) the Polygon is invalid GeoJSON.

I am looking for a way to integrate cleanup/fixing functionality for my user-generated Polygons in my Node.js application that uses a MongoDB (Turf being my main tool for working with GeoJSON).

I have briefly looked at ST_MakeValid and prepair which both seem to do what I need, but I don't see an easy way to integrate either into my stack.

The question is whether there are any libraries that are either written in JavaScript or can easily be wrapped to be used via an npm and play nicely with my existing stack and doesn't require bringing in heavy-weights like PostGIS.


texts:
What are the uses (public and private industry) of building footprints?

I'm trying to figure out how polygon building footprints can be used.  i.e. are they used for determining tax liabilities, etc...


texts: Clipping vector layers with different CRS

texts:
I have two vector layers:

provinces (Belge_Lambert_1972, EPSG:103300) - clip layer
building (WGS 84/ Pseudo Mercator, EPSG:3857) - input layer

I would like to overlay them, however there is something wrong. I believe it had to do with the Coordinate Reference System.

When i work with the layer without conversion of the reference systems, the data does overlay exactly. I know that the buildings should lie within the provinces from the definition of the layer. However, when I want to clip the layers I get an notification:

  Layers do not all use the same CRS. This can cause unexpected results.
  Do you want to continue?

On the other hand, when I do convert the WGS 84/ Pseudo Mercator, EPSG:3857 to Belge_Lambert_1972, EPSG:103300 CRS, the layers are not overlapping at all. Actually, I can find the layers by using "zoom to layers".

How can I clip these two layers correctly?


texts: Create metadata for multiple datasets simultaneously in ArcGIS 10

texts:
I am trying to find a way to create metadata through Arc Catalog for multiple datasets simultaneously. I'm working within a file geodatabase, which composes of multiple feature classes and rasters. The data has common theme so I'd like to create template of metadata for one item and then populate other datasets within the same metadata information. I'm planning to use ISO standard metadata template.

I've done some research and it appears that there are three options (but neither of them too successful so far):
 - use 'Metadata Importer' tool within Conversion Toolbox (however I will be able to copy metadata on one-to-one basis only)
 - X-tools pro have 'Batch Metadata Editing' option available but after setting the parameters within ArcMap I can't see these canges through Arc Catalog
 - third option is adding 'Batch Importer' (http://edndoc.esri.com/arcobjects/9.0/Samples/Metadata/Importers/Batch_Importer/Batch_Importer.htm). I managed to add it to ArcCatalog but I'm not able to run it (I get error message 'ICommand_OnClick() ).

I'm working with ArcGIS 10.


texts: QGIS: Turning a string into "relative" hyperlink

texts:
I have a few columns in my shapefile that to refer to certain folders and files (PDF, JPG etc.). Because these are interesting in a form for each object, so that people can open them. I want to make a "relative" hyperlink to those files, and by relative I mean relative to where the shapefile itself is.
For example:
xxx/Shapefile folder
xxx/shapefile folder/Folder 1/File 1
xxx/shapefile folder/Folder 1/File 2
xxx/shapefile folder/Folder 1/File 2

In the atrribute table is this: /Folder 1/File 2
How do I do this?


texts: QGIS + PostgreSQL - file that works on one PC will not work on another

texts:
I have QGIS on 2 different standalone PC and I share/sync files between them.
I now need my second PC QGIS to work with a file that includes PostgreSQL queries.
PC 1 is original PC and all is well there.  Various versions of QGIS and all deployed versions of postgres/postgis behave as expected.  All hosted on standalone PC (localhost)
PC 2 - QGIS works fine, I have created an identical Postgresql database including GIS extensions, name, credentials etc and populated this with identical underlying data.  QGIS can connect and query this database.

Issue: If I then open one of my sync'd QGIS files on PC 2 it fails to connect/open any vectors linked to the PostgreSQL database.  I've tried to troubleshoot/google, used various V3 QGIS and failed to resolve this. Same file opens fine on PC 1. On the face of it I can see no reason this won't work - however need some expert opinion/advice.


texts: Raster clipper is giving negative coordinates in QGIS?

texts:
I have an ECW file that I want to clip and to convert to MBTILES. It's on ITM CRS, and when I use clipper on QGIS the output is in the same coordinates but negative! So 250000, 750000 becomes -250000,-750000

Here's the printscreen


texts: Same Origin Policy Problem using WFS-T with Geoserver & OpenLayers 3

texts:
I'm trying to use a WFS-T Service to insert a simple Point from OpenLayers 3 into Geoserver in a VirtualBox (test environment). There is always the same error message:

    Failed to load resource: Origin http://localhost:8383 is not allowed by Access-Control-Allow-Origin. at http://HOST_IP:8082/geoserver/test/wfs
XMLHttpRequest cannot load http://HOST_IP:8082/geoserver/test/wfs. Origin http://localhost:8383 is not allowed by Access-Control-Allow-Origin. (13:29:51:925 | error, javascript) at resources/webappv1.html

I took the issue and tried to solve it with Cross-Origin Resource Sharing (CORS).
(http://en.wikipedia.org/wiki/Cross-origin_resource_sharing)

A ProxyPass (like http://www.gistutor.com/geoserver/21-intermediate-geoserver-tutorials/38-configuring-geoserver-proxy-for-public-and-remote-data-access.html)  would be for security reason only a bad alternative.

I found a solution on the Geoserver Mailing Archive by Dusan Fedorcak (Search for [Geoserver-users] CORS for jetty 6.1.8 (Geoserver 2.x)) and tried to manage it , but it is still not working.

Where do I go wrong?
Or is there an other workflow to make a post Request with OL3?

Here my versions:

OpenLayers 3
Geosever 2.4.4 runs on VirtualBox 4.3.10 with OSGeo Live7.9 installation
Jetty 6.1.8
PostGIS 2.1.1

Here my code:

var xmlData =
            +'< wfs:Transaction\n'
            + 'service = "WFS"\n'
            + 'version = "1.0.0"\n'
            + 'xmlns:wfs = "http://www.opengis.net/wfs"\n'
            + 'xmlns:gml = "http://www.opengis.net/gml"\n'
            + 'xmlns:test = "http://test.org"\n'
            + 'xmlns:xsi = "http://www.w3.org/2001/XMLSchema-instance" \n'
            + 'xsi:schemaLocation = "http://www.opengis.net/wfs \n'
            + 'http://HOST_IP:8082/geoserver/schemas/wfs/1.0.0/WFS-basic.xsd \n'
            + 'http://test.org/ \n'
            + 'http://HOST_IP:8082/geoserver/test/wfs/DescribeFeatureType?

typename=test:user_poi">\n'
            + '< wfs:Insert >\n'
            + '< test:user_poi >\n'
            + '< test:gid >' + gid + '< /test:gid>\n'
            + '< test:id >' + id + '< /test:id>\n'
            + '< test:tid >' + id + '< /test:tid>\n'
            + '< test:auid >' + auid + '< /test:auid>\n'
            + '< test:mdid >' + mdid + '< /test:mdid>\n'
            + '< test:text >' + text + '< /test:text>\n'
            + '< test:the_geom >\n'
            + '< gml:Point srsDimension = "2" srsName = "urn:x-ogc:def:crs:EPSG:3857" \n'
            + '< gml:coordinates decimal = "." cs = "," ts = " " >' + lat + ',' + long + '< /gml:coordinates>\n'
            + '< /gml:Point>\n'
            + '< /test:the_geom>\n'
            + '< /test:user_poi>\n'
            + '< /wfs:Insert>\n'
            + '< /wfs:Transaction>\n';

    $.ajax({
        type: "POST",
        crossDomain: true,
        url: wfsurl,
        dataType: "xml",
        contentType: "text/xml",
        data: xmlData,
        xhrFields: {
            withCredentials: true
        },
        username: 'user',
        password: 'user',
        success: function(data, textStatus, request) {
            alert(request.getResponseHeader('some_header'));
        },
        error: function(request, textStatus, errorThrown) {
            alert(request.getResponseHeader('some_header'));
        }
    });

Thanks for your help!


texts: Dividing up large region into tiles and visualizing them (Descartes Labs)

texts:
I would like to run a deep learning model over an AOI that is about 300,000 square kilometers. I used the dl.scenes.DLTile.from_shape function to automatically divide the region into tiles.
Is there an easy way to visualize the grid of tiles that was generated?
Just the polygons themselves, not having to pull any imagery.


texts:
I have some OGC WMS services made in MapServer. Now, I want to restrict access to data for some users by, for example, username, number of requests, or bbox.
Can this be done using MapProxy?
I read the documentation and there were mentioned WSGI middleware and MapProxy Authorization API as ways to accomplish this, but I am having trouble implementing this.
I am using MapProxy inside a docker container.


texts: Leaflet geoJson mouseover does not work for all polygons

texts: I have  geojson polygons. Code highlights the polygon on mouse over.  The problem is: for some polygons, the mouse over does not activate; and sometimes I need to move mouse outside polygon and back in a few times to get it to activate.


texts:
I am having difficulty installing QGIS on my Mac OS X 10.7.5 (Lion). I understand that QGIS requires GDAL Complete 1.10 framework package and Matplotlib Python module. I installed GDAL (which includes PROJ, UnixImageIO, GEOS and SQLite3) and then I installed NumPy.

When I try to install QGIS, a message pops up saying that:

  QGIS requires the Matplotlib python module

I thought I took care of that by installing NumPy.

Does anyone have any advice to share?


texts: Using CARTOframes to push to CARTO account

texts:
I am using CARTOframes to create and style a map that includes widgets. I'm able to push a Dataset to my Carto account and to publish and update a map, however I can't seem to get that published map to appear in my dashboard when I login to Cartodb. The goal is for people on my team to be able to make changes to a map either from a jupyter notebook using something like CARTOframes OR through logging on and using the Builder manually.

I have version 1.0b2 of CARTOframes installed and here is the code:

#setting credentials:
set_default_credentials('https://{}.carto.com/'.format(username), api_key)

Creds = Credentials(username=username, api_key=api_key)

#to turn geopandas geodataframe into a carto Dataset:

locations_data = Dataset(df)

#to create map

testmap = Map(Layer(locations_data,
          popup=Popup({'hover':[{'title':'service', "value":'$service_name'}]}),
            widgets=[
            default_widget(
                title='NYC Social Services Locations',
                description='Visualization of service sites for city-funded social services.',
                footer='Data source: SSSLD'
            ),
            formula_widget(
                'count',
                is_global=True,
                footer='sites'
            ),
            category_widget(
                'org_name',
                title='Provider Organization Name',
                description='Select an org to filter',
            )
            ]
         ) )

#to sync data with my account:

testmap.sync_data('test_data')

#to publish map

testmap.publish('sssld_test_map_jupyter', maps_api_key=api_key, credentials=Creds)

#to update published map:

testmap.update_publication('sssld_test_map_gui', password=None, maps_api_key=api_key)


texts: QgsVectorLayer creation in Qthread (pyQt)

texts:
Our QGIS project is dynamically generated on startup based on user roles. There are nearly 120 QgsVectorLayers that need to be generated. If created sequentially, the login process takes up to 2 minutes which is really annoying.
So I thought I'll use a QThreadPool and create the layers in there. So the threads basically just call this kind of code:
class LayerLoaderWorker(QRunnable):

    def run(self):
        db_connect = DatabaseConnection()
        uri = QgsDataSourceUri()
        uri.setConnection(db_connect.host, str(db_connect.port), db_connect.dbname, db_connect.user,
                          db_connect.password)
        uri.setDataSource('public', self.viewname, 'geometry', '', 'id')

        locker = QMutexLocker(mutex)
        self.layer = QgsVectorLayer(uri.uri(False), self.layerSettings.displayname, "postgres")
        self.finished = True
        self.loader.workerFinished()
        locker.unlock()

Only the layers are created in the threads, after creation they are added to the QgsProject in the main (UI) thread.
Sometimes it seems to work, but most of the time the layers are not loaded correctly and the geometry column is empty.
Is the creation of the QgsVectorLayer object thread safe? It seems to me that this isn't the case. I read somewhere that UI-related things should be handled in the UI-Thread, but I thought creating the object (which seems to connect to the database which takes some time) is a non-ui-task.
Is it possible to load layers in parallel somehow?


texts: Has anybody used arcpy.mapping.ExportReport successfully on 64bit machine?

texts:
The online help for arcpy.mapping.ExportReport says that it has some restrictions on use but does not list 64bit machines amongst them.

I am running the following test from IDLE under ArcGIS 10.1 SP1 (with 64bit Background Geoprocessing installed) on Windows 7 SP1.

I've create a simple test.mxd (one layer called Suburbs in data frame called Layers) and a simple test.rlf (included Suburbs, said Finish to take all defaults, saw report and saved it as *.rlf)

However, this code (more or less straight from the Help example):

import arcpy
mxd = arcpy.mapping.MapDocument(r"C:\temp\test.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
lyr = arcpy.mapping.ListLayers(mxd, "Suburbs", df)[0]
arcpy.mapping.ExportReport(lyr,
                           r"C:\temp\test.rlf",
                           r"C:\temp\test.pdf")
del mxd

throws this error:

Python 2.7.2 (default, Jun 12 2011, 14:24:46) [MSC v.1500 64 bit (AMD64)] on win32
Type "copyright", "credits" or "license()" for more information.
>>> ================================ RESTART ================================
>>>

Traceback (most recent call last):
  File "C:\temp\test.py", line 7, in <module>
    r"C:\temp\test.pdf")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\utils.py", line 181, in fn_
    return fn(*args, **kw)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\mapping.py", line 515, in ExportReport
    return report_source._arc_object.ExportReport(*gp_fixargs((report_layout_file, output_file, dataset_option, report_title, starting_page_number, page_range, report_definition_query, extent, field_map), True))
RuntimeError: Error in generating report
>>>

Any thoughts on cause and/or possible workaround are welcome!
(My suspicion is that ExportReport may not be supported on 64bit machines.)


texts:
I'm classifying images of Red River Delta. Unfortunately, the region is split into two consecutive Sentinel 1A images. Therefore, I want to mosaic them into single image for convenient classification.

How do I do this?

Here is the code:

// Use these bands for prediction.
var bands = ['VV', 'VH'];

// Load a sentinel 1A image to be used for prediction.
var img1 = ee.Image('COPERNICUS/S1_GRD/S1A_IW_GRDH_1SDV_20170104T225034_20170104T225059_014688_017E5C_6A39').select(bands);

var img2 = ee.Image('COPERNICUS/S1_GRD/S1A_IW_GRDH_1SDV_20170104T225059_20170104T225124_014688_017E5C_750B').select(bands);

var img_temp = ee.Image.cat([img1, img2]);


texts: Do I need to use ST_SetSRID, or does ST_MakePoint()::geography work?

texts:
According to this SO answer, the way to convert a lat/long to a geography type is

ST_SetSRID(ST_MakePoint(long, lat), 4326)::geography

but so far I've just been using ST_MakePoint(long, lat)::geography and it seems to work fine. My queries include checks for

ST_Covers(shape, ST_MakePoint(long, lat)::geography)

and I haven't seen that give an incorrect answer. But that doesn't mean it's always going to work.

The shape column is defined as

shape geography(MULTIPOLYGON, 4326)

so I thought maybe having the SRID in the column definition was doing the same work as ST_SetSRID. But I really don't know.

Probably not relevant, but just in case: the data comes from Natural Earth's countries data, converted to SQL with

shp2pgsql -GD ne_10m_admin_0_countries.shp table_name

That creates a table with a geog column also defined as geography(MULTIPOLYGON, 4326), and I've copied some of the data from that into the table with shape.


texts: Incorporating custom/image markers into Leaflet MarkerCluster?

texts: I'm using Leaflet MarkerCluster. It works fine, however, when you drill down to an individual marker (not cluster) I can't figure out how to replace the standard marker with my own custom icon/image.


texts:
I'm a software developer, don't know much about the hardware specifications and don't know where to start. (Even not sure is that the right stackexchange I should ask this question.)

So, I want to know about (1) the cost and (2) document, link, etc. for the following specifications:

I want to make a project which I will put a chip (?) to a user and track that users 3D location with at least 10 cm accuracy, 10 times a sec, in a particular approximately 25x40 meters place and use that 3D location in my C, C++, Python, Java application. I read something about triangulation, and it makes sense. However that's all I know. I don't know any device that I can buy to do such...

Can you send me something to read?


texts: Finding Network Data for Route Optimisation in Saudi Arabia using ArcGIS?

texts:
I have to do a route optimisation through ArcGIS vrp solver for 1200 pizza drop points/daily, for a city based in Saudia Arabia.

Where can i get the road  network dataset?


texts: Calculate Blue sky albedo from Black sky albedo and White sky albedo

texts:
I have got values for Black sky albedo and white sky albedo, from MODIS imagery. What is the simplest way to calculate blue sky albedo? I want to use albedo for snow cover.

A website which has listed this formula:

Blue-sky (Actual) Albedo = White-sky Albedo * f(optical depth, solar zenith angle, aerosol type, band) + Black-sky Albedo * (1 - f(optical depth, solar zenith angle, aerosol type, band))

Now the thing is I don't have solar zenith angle, aerosol type etc..

What do I do? Please help


texts: How to use nested for loop to iterate over features of a layer?

texts:
I have a shapefile having features with duplicate IDs. I want to merge features with same IDs and take average of values of other field which are different. How is it possible in python in QGIS3?

Here, all other fields have same values of both features except last two fields. I want to merge the features and take average of values of last two fields. Here ID is same but last three fields are different because the file is the output of intersection, where one feature intersects with more than two features of other file so it is joining attribute of overlayed shapefile but input shapefile is repeating the ID and dividing a feature into two different features.
I have tried the following code but it is not working.

for f in interLayer.getFeatures():
        g =next(interLayer.getFeatures())
        if f['PIPE_ID']!=g['PIPE_ID']:
            f['k-new']=f['k']
            interLayer.updateFeature(f)
        if f['PIPE_ID']==g['PIPE_ID']:
            f['k-new']=(f['k']+g['k'])/2
            interLayer.deleteFeature(g.id())
            interLayer.updateFeature(f)
    interLayer.commitChanges()

Instead of delete, I want to merge features. But the code is not working for delete too. Here k-new is the field created to append averaged value.also, the duplicate ID can be more than two in some cases.


texts: Setting vertical scan angles in visibility analysis (viewshed) QGIS?

texts:
Is it possible doing a visibility analysis (viewshed) in QGIS setting a vertical angle of search from the observer points?
In ArcGIS Pro viewshed, there are the vertical upper angle and vertical lower angle parameters. In Global Mapper viewshed, it is possible to set angles in the tool dialog, 'Receiver Elevation' section.
In QGIS, I could not find such option in the Visibility Analysis plugin (v. 1.2) and in the r.viewshed tool from GRASS (gdal_viewshed does not have it either).
Is there a tool or a way I could work around on this, while getting same results?
PS: I am aware about the target_elevation parameter (which is the one I am currently using), but it is a different (substitute or complementary) method from setting a vertical angle range from the observer view.


texts:
Running ArcMap 10.1

So, I have a bunch of layers that all have the same name. I'm trying to rename them so that they have names that actually mean something.

Everything runs fine until it is time for the names to actually change. Nothing happens, and I don't know why.

Here is my code:

import arcpy
mxd = arcpy.mapping.MapDocument(r"C:\Users\User\Documents\GISProj\Map6.mxd")
names=['Alexander','Avery','Buncombe','Burke','Caldwell','Cleveland','Durham','Gaston','Henderson','Jackson','Macon','Orange','Polk','Rutherford','Stokes','Surry','Transylvania','Watauga','Wilkes']
names.reverse()
cntys=arcpy.mapping.ListLayers(mxd,'slope')
for c in cntys:
    c.name=names[cntys.index(c)]
arcpy.RefreshTOC()

after I run it, I can see that cntys returns

[<map layer u'Wilkes'>, <map layer u'Watauga'>, <map layer u'Transylvania'>, <map layer u'Surry'>, <map layer u'Stokes'>, <map layer u'Rutherford'>, <map layer u'Polk'>, <map layer u'Orange'>, <map layer u'Macon'>, <map layer u'Jackson'>, <map layer u'Henderson'>, <map layer u'Gaston'>, <map layer u'Durham'>, <map layer u'Cleveland'>, <map layer u'Caldwell'>, <map layer u'Burke'>, <map layer u'Buncombe'>, <map layer u'Avery'>, <map layer u'Alexander'>]

just like it should. But the actual layer names in the TOC don't change, even when I save, close, and reopen the map.

From what I'm looking at here:https://geonet.esri.com/thread/43148
and here:Use ArcPy to rename MXD layers
it should work, but it doesn't. Whats up with this?


texts: Correctly displaying text in QGIS that is sourced from a PostGIS database

texts:
Long time listener, first time caller...

I have an issue with displaying the text that is called from a PostGIS 2.4 database into QGIS 2.16.1.

The text is displaying in the table like so; {ACACIA AVENUE} but when it displays in QGIS it does so like; {"ACACIA AVENUE"}.

To resolve this I have labelled this column using the following in the 'Label with' box in Layer Properties using substr:

CASE WHEN  "streetType" = 'Numbered Street' THEN concat("nationalRoadCode", ' ', substr( "name",3,((length("name")-4)))) ELSE substr( "name",3,((length("name")-4))) END

Which works absolutely fine and ACACIA AVENUE displays as it should!!

However, I subsequently find that where there is a single street name in this particular column it is called across into QGIS like this; {MERESIDE} so when this is labelled using the code above it comes out as ERESID!

I would rather resolve this issue using styling as opposed to altering the table as this would mean something else to do when the table is refreshed with CoU.

So my question is, how can I display both the double word texts and single word texts without losing the leading and trailing characters?


texts: How to index mbtiles to speed up load time?

texts:
I have created a mbtile based on OpenStreetMap data with Maperitive using the command:

generate-mbtiles file=Path\filename.mbtiles maxzoom=16

Loading the file in QGIS 2.14.0 and zooming takes 15-30 seconds, even if zoomed in at zoomlevel 16.

Is there a method index mbtiles to speed up load time?


texts: How do I use "execute SQL" from QGIS graphical modeler?

texts:
I have a points layer representing road gullies and a non-spatial table of gully inspection data, both stored in a postGIS database. I would like to create a QGIS model to perform some repetitive analysis tasks on the tables, once they are joined on a common id field.

I have tried to use the "join attributes table" tool in the model, but there seems to be limited control over how the join is performed. For one of the analysis tasks I need just the most recent inspection, where there may be 3 or 4 inspections linked to the same gully.

I then found the "execute SQL" tool which can be called by the model. However, I've tried multiple versions of syntax and I just can't seem to get this tool to work for me. Even when executing the tool from outside the model using just

SELECT * FROM gullies

I get an error saying "Unable to open datasource ... with the following drivers..."

Can someone please let me know how I could perform this type of join in the graphical modeler?


texts:
I have 4 shapefiles of India state from 4 decades, and I trying to identify changes within the sub-district level such as splits or merge of sub-dist, border changes etc.

I am using R but I can use other packages such as ArcGIS...

Any suggestion how to check and measure the border change over time?


texts: Table to dBase tool failure when converting .xls files

texts:
Whenever I load my Excel 2010 worksheets into ArcMap 10 Table to dBASE tool, it fails, returning the error: "Failed to Copy Rows" The Copy Rows tool works just fine for me on a single document. There are no special characters or spaces in my file paths, workbook or worksheet names, or in column headers.

Has anybody else out there had an issue with this tool? Or is there a simple fix for this that I am missing?


texts: Editing shp layers style simultaneously in QGIS?

texts:
I have some polygons in a shp file that I need to change their line width.
I head joined them in one group

Is there a way to change the line with for all of the layer in the group simultaneously and not one after the other?


texts: Mask out areas inside a feature collection on Google Earth Engine

texts:
I am quite new to javascript and I am trying to apply a mask on an image to exclude areas inside a feature collection, I have seen that for the case of image collections the "mask" and "mask.not" expressions have opposite effects, however for "clip" that is not the case, then what would be the opposite expression to clip?

Here is my code with the mentioned mask missing, I would like to mask out the protected areas:

///load collections country boundaries
var country =ee.FeatureCollection('USDOS/LSIB/2013');

// load protected areas
var PA= ee.FeatureCollection('WCMC/WDPA/current/polygons');

// load elevation and mask for a country
var srtm = ee.Image('USGS/SRTMGL1_003').clip(country.filter(ee.Filter.eq('cc', 'CO')));
Map.addLayer(srtmmask, {min: 1, max: 100, palette:['green', 'yellow', 'red']}, 'Elevation');


texts: Comparing differences between calculated areas of multiple shapefile by using PyQGIS and geopandas python module

texts:
I have calculated area of multiple shapefile by using PyQGIS and geopandas python module but results are not same. Which one is better?

logic in QGIS:

with edit(eligibleLayer):
    for gFeat in eligibleLayer.getFeatures():

        area = QgsDistanceArea()
        area.setEllipsoid('WGS84')
        area.setEllipsoidalMode(True)
        area.computeAreaInit()

        geom = gFeat.geometry()
        print("Feature ID %d: " % gFeat.id())
        attrs = gFeat.attributes()

        multiPart = gFeat.geometry().isMultipart()

        if (multiPart is False):
            pol = gFeat.geometry().asPolygon()
            print("glac_name"+attrs[14]+'-'+"area is %.2f km2" % (area.measurePolygon(pol[0])/1e6))
        else:
            pol = gFeat.geometry().asMultiPolygon()
            print("There are %d polygons" % len(pol))
            i=1
            for geom in pol:
                print("area%d is %.2f km2") % (i, area.measurePolygon(geom[0])/1e6)
                i +=1

GeoPandas with python:

fp = gpd.read_file('/home/priya/project/EAST/data/east_him.shp')

a = fp.crs
tost = fp.copy()

tost = tost.to_crs({'init': 'epsg:3857'})
b = tost.crs
print(b)

tost["su_area"] = tost['geometry'].area / 10**6
tost.head()


texts: Why is "Selected Graphics (Clipping)" always grayed-out in Export Raster Data?

texts:
I would like to utilize the "Selected Graphics (Clipping) option so I may export my raster data along the boundaries of a polygon mask which I created earlier.

I would use "Selected Graphics (Clipping) to do this, however it is always grayed out.

Why is this?


texts:
I'm trying to populate a field (AVO) using this statement

AVO =
dim avo
if [join] = 1 then AVO = "YES"

It works if I do it in ArcMap but it doesn't when I use the same syntax in ModelBuilder.


texts: How to only check one layer in a group at a time?

texts:
I have a group of layers and in my plugin, I would only like one layer from a specific group to be checked at any time. The plugin contains a toggle button which connects each layer to a function. The function finds the checked layer and adds it to a list, it then finds all other layers and I tried to set their node state to Unchecked but clearly I'm missing something.

Here is the code:

root = QgsProject.instance().layerTreeRoot()
group = root.findGroup('Main group')

def check():
    if self.plugin_button.isChecked():
        # Create empty list to contain checked layer
        layer_list = []
        for child in group.children():
            node = root.findLayer(child.layer().id())
            if node.isVisible():
                # Add checked layer into list
                layer_list.append(node.name())
            # For layers not in list
            if child.name() not in layer_list:
                ###### Not sure what to do here ######
                node = root.findLayer(child.layer().id())
                node.setVisible(Qt.Unchecked)
                iface.legendInterface().setLayerVisible(child.layer(), False)
    else:
        for layers in group.children():
            # Disconnect layers in group from function
            layers.visibilityChanged.disconnect(check)

# Connect layers in group to function
for layers in group.children():
    layers.visibilityChanged.connect(check)

Update:

Made some progress, the following code now only allows one layer to be active at any time. Only issue is when a layer is checked, clicking on any layer unchecks all layers. Still looking to see how to keep one layer checked at all times. Also if a layer is checked, I don't want it unchecked (i.e. I want it 'locked'). The only way to uncheck it is to select another layer.

root = QgsProject.instance().layerTreeRoot()
group = root.findGroup('Main group')

def check(layer):
    if self.plugin_button.isChecked():
        layer_list = []
        layer_list.append(layer.name())
        for child in group.children():
            if child.name() not in layer_list:
                iface.legendInterface().setLayerVisible(child.layer(), False)
    else:
        for layers in group.children():
            layers.visibilityChanged.disconnect(check)

# Connect layers in group to function
for layers in group.children():
    layers.visibilityChanged.connect(check)


texts: Partial loss of image when exporting from composer in QGIS

texts:
When I Have a map in composer this is what I see

but when I export as an image or PDF this is what I get

Wondering why a small part of the image is missing on the left side of the page and how I can correct this?


texts: Identify all layers on the map without pre-designing the infoWindow

texts:
This applies to the ArcGIS Server JavaScript API 2.4.

The Identify Task sample shows how IdentifyTask can be used to find features in multiple layers, as specified in the IdentifyParameters:

identifyParams.layerIds = [0,2];

In this sample, the infoWindow already contains  elements for the expected layers, which are used to display multiple features in separate tabs of the infoWindow:

//In the JavaScript:
dijit.byId("bldgTab").setContent(layerTabContent(bldgResults,"bldgResults"));
dijit.byId("parcelTab").setContent(layerTabContent(parcelResults,"parcelResults"));

Is there a way to handle this dynamically, and create the required tabs at run-time?

For example, I wish to run IdentifyTask on all visible layers, and return the results in a separate tab of the infoWindow. How can I create the tabs dynamically?


texts: FME workflow not outputting tab file correct

texts:
My FME workflow, combines property data sourced from three different datasets.  It must start and end with the same number of features, so I have in place a statisticsCalculator at the start and finish of the workflow, added a tester to ensure they are the same, and added another tester to ensure the count is above 24,000 (which I know the output will never be less).

Now very sporadically neither of these safety nets work and the layer is still produced.  Usually we can fix the issue in 10min the next day but its frustrating I can't figure this out or why it happens sporadically, and its hard to test any improvements.

I have reviewed the log file, and get these errors when I have an issue.  The file gets overwritten every night to catch any updates to databases so there is adequate disk space, the file isn't in use by another application or workflow, and the errors regarding attributes - the file is written fine majority of the time.

Error 1

2016-07-24 19:27:45|  92.1|  0.0|ERROR |MapInfo file 'C:\etc' could not be created.  Check that file path name is valid and there is adequate disk space. Also ensure that the total record length for all the attributes is less than 4000 bytes and number of attributes are less than 250

Error 2

MapInfo file 'C:\etc' could not be created.  Check that file path name is valid and there is adequate disk space.

Error 3

Also ensure that the total record length for all the attributes is less than 4000 bytes

With the implemented precautionary steps not capable of capturing the failure, is there another way to add a safety net actually into the writer to ensure that it is not written if say it doesn't have more than 24,000.  The image below shows the workflow.


texts: Why does GeoNames say there are only 73 parks in Mexico?

texts:
I did this GeoNames query

http://api.geonames.org/search?featureCode=PRK&username=genadinik&country=MX&style=full&maxRows=10

and it says there are only 73 parks in all of Mexico. Am I doing something wrong, or GeoNames mostly covers the United States?

Thanks!


texts: I have a QHD monitor and it's been giving me all sorts of problems when I build my custom plugins in QGIS3. I've built my plugins using the plugin builder and designed the UI with QT Designer for QGIS - but on monitors at certain resolutions the texts get bunched up and nothing I do in the QT Designer interface fixes it. Is this potentially a problem on the QGIS-side or solvable in QT?


texts:
I have developed a QGIS plugin that is useful for many people and it reduces a day-long workflow to a few minutes. I have ideas for other plugins too. To continue to add value to QGIS and save the time of other users, I would like to know if I can earn money through QGIS plugins.
Can we put pop-ups in the plugin that every now and then ask for user feedback and if they like the plugin then ask them for donation or is there a better way than this to raise money through the QGIS plugins?


texts: Multiple date filter Google Earth Engine

texts:
I am trying to filter an ImageCollection by multiple time intervals (rain seasons) that change from year to year. Hence, I can't just use the ee.Filter.calendarRange.

I tried to merge the dates together and filter the collection but it does not work. What do I have to do?

var ROI = /* color: #0b4a8b */ee.Geometry.Polygon([
      [36.5048217773437, 0.560293804172086],
      [37.9302978515625, 0.560293804172086],
      [37.9302978515625, 1.875580463973453],
      [36.5048217773437, 1.875580463973453]
  ]);
    Map.centerObject(ROI,8);

var s2 = ee.ImageCollection('COPERNICUS/S2')
  .filterBounds(ROI);

var longRain2016 = s2.filterDate('2016-04-01','2016-05-31');
var shortRain2016 = s2.filterDate('2016-11-15','2016-12-31');
var longRain2017 = s2.filterDate('2017-04-01','2017-05-31');
var shortRain2017 = s2.filterDate('2017-10-01','2017-12-31');
var longRain2018 = s2.filterDate('2018-03-01','2018-05-31');
var dateFilter = ee.ImageCollection(longRain2016.merge(shortRain2016).merge(longRain2017).merge(shortRain2017).merge(longRain2018));

var s2dateFilter = s2.filter(dateFilter);

var count = s2dateFilter.size();
print('Number of filtered images in the region of interest: ', count);


texts: CentOS 6.5 - After compiling gdal.so: Unable to load dynamic library

texts:
CentOS release 6.5 (Final)
Linux version 2.6.18-028stab101.1 (root@rhel5-build-x64) (gcc version 4.1.2 20080704 (Red Hat 4.1.2-46)

I'm trying to Install Cartaro CMS (http://cartaro.org/). GeoServer, Postgresql and PostGIS are up and running.

The only missing software requirement to install Cantaro is the GDAL-PHP extension. I've spend hours to resolve the dependencies but after successfully compiling the gdal.so, "make test" fails with the following warning:

PHP Warning:  PHP Startup: Unable to load dynamic library '/tmp/php5-gdal-0.1.3/modules/gdal.so' - /tmp/php5-gdal-0.1.3/modules/gdal.so: undefined symbol: _ZN17GDALDriverManager10MoveDriverEP10GDALDriveri in Unknown on line 0

PHP_SAPI    : cli
PHP_VERSION : 5.4.36
ZEND_VERSION: 2.4.0
PHP_OS      : Linux - Linux 268.xxx.de 2.6.18-028stab101.1
INI actual  : /tmp/php5-gdal-0.1.3/tmp-php.ini
More .INIs  :
CWD         : /tmp/php5-gdal-0.1.3
Extra dirs  :
VALGRIND    : Not used
SKIP Test: OGRGeometry [tests/ogrgeometry.phpt]
FAIL Test: Version number [tests/version.phpt]

I've compiled the actual stable vesions of gdal (1.9.2) and proj4 (4.8.0) by following these instructions:

https://github.com/wavded/ogre/wiki/Compiling-a-recent-ogr2ogr-from-source-on-CentOS-%28RHEL%29


texts: Calculate centroids in memory from polygon using PyQGIS

texts:
I need to calculate Centroids from polygon layer without reading and writing shapefiles.
I tried this code in python console (QGIS 3):

layer = iface.activeLayer()
params = {'INPUT': QgsProcessingFeatureSourceDefinition(layer.id(), True),'ALL_PARTS':False,'OUTPUT': 'memory:'}
newLayer = processing.run("native:centroids", params)
layerOutput = newLayer['OUTPUT']
QgsProject.instance().addMapLayer(layerOutput)

There are no error's messages and the new layer has the right fields, but there are not features, it's empty. If I run the algorithm using processing toolbox it works.
Where is the error?


texts: Allow client to upload ESRI shapefile in OpenLayers 3

texts:
I want to build functionality to upload ESRI shapefile and overlay it on map using OpenLayers 3. Once user cross checked I want to save to PostgreSQL/PostGIS table.

So far as I Googled I am not getting any resources. Is there any way to do in OpenLayers 3?


texts: Leaflet Routing Machine routing icons are not displaying correctly!

texts:
I'm having an issue with the routing icons in Leaflet Routing Machine. For some reason that I can't figure out, the only icon that is displayed is the up arrow icon, regardless if the router says to turn left, right, or etc. I've followed the specifications and still am having this issue.

Any clues as to what could be the cause?

Here are some screenshots to display my issue:


texts:
I haven't been able to successfully use the QGIS MetaTools plugin for describing metadata. Has anyone used it before?
When I try to create metadata for a shapefile, the plugin creates the xml file but then I get a python error:

An error has occured while executing Python code:

Traceback (most recent call last):
File "C:/Users/daniel/.qgis2/python/plugins\metatools\metatools.py", line 267, in doEdit
standard = MetaInfoStandard.tryDetermineStandard(self.metaProvider)
File "C:/Users/daniel/.qgis2/python/plugins\metatools\standard.py", line 35, in tryDetermineStandard
print "TryDetStd ", text, " Text"
IOError: [Errno 9] Bad file descriptor

I'm using QGIS 2.4.0-Chugiak Chugiak, installed from the OSGeo4W installer with python version 2.7.4. Has anyone seen this error?


texts: FME: Connecting the output of a workflow to the input of another workflow

texts:
Solved: I've tried the advice from @Mapperz - Connecting the transformer before the Writer is writing the shp with the transformers "AttributeVreator_2 & Dissolver".
  Anyway, I'm still interested on my question. It is possible to connect a Writer and Reader module, that the Reader uses the output of another workflow?

Version: FME Desktop 2014 SP3

I have 2 separate workflows which I run one after other. This is a litte bit inefficient, therefore i want to optimize the process by  connecting the two workflows.

Workflow 1 generates a shape with lots of attributes as output.
Workflow 2 uses this output (shape) to do more operations.

Problem: I'm unable to link the output of Workflow 1 to the input of Workflow 2. I tried it with the Transformer "WorkspaceRunner", but it is not possible to define the right input for me.

Maybe there is antoher way to solve this problem?


texts: How to extract nodes with PostGIS to a multipoint geometry

texts:
I have a line in a PostGIS table test_line. I want to extract it's nodes into a multipoint geometry.

I tried the following:

SELECT ST_Collect((ST_DumpPoints(the_geom)).geom) FROM test_line;

But it gives me the error ERROR:  set-valued function called in context that cannot accept a set.

What do I have to change in order to get it to work?


texts: Embedding qgis2web exports into WordPress?

texts:
I exported my QGIS project with the qgis2web plugin as OL3.

Can I embed this folder to a wordpress homepage or can I import it to google maps and add this to WordPress?


texts: How to join several inner and outer lines that share the same field in a single polygon?

texts:
I have a data source with thousands of lakes, ponds and small puddles.

Using an algorithm of anitagraser.com I have managed to clean the ponds and small puddles, but with the big lakes I have a problem. They are formed by many lines, and none runs completely by its edges, are only pieces.

I put an image where you can see only one of this lakes in full:

And another image where you can see (in different colors) some of the lines that form it (they are not all, only a few, too much work mark each line in a different color, but you can have an idea)

I've tried playing with 'dissolve' and with v.clean.snap, but I can not make much progress.

Do you have any idea how I can manage to make all this tangle of lines 1 only polygon? It must be a process that can be automated, nothing manual, since there are thousands of lakes like this one into the source. I am able to make complicated scripts, but I am an inexperienced one in GIS (using QGIS and PyQGIS). (Added after Evan comment -> The details are not important, it is more important to close the polygons, I would sacrifice the details if you can close the polygon and create a single figure)

Yes, we have something on our side, all the lines that form a large lake have the same label in the attribute table.

You can download this shape here: big_lake.zip

Note: My project (map) will be a open project, and all data will be save into OSM if possible.


texts:
I have got a shapefile, and a GeoTIFF image. I want to crop my GeoTIFF based on my shapefile. EASY ! but I don't know I get lots of error.

I tried

library(raster)
xp <- raster("c:/test/M145.tif")  ## 1
shp <- shapefile("c:/shapefile/Indy_Bdry.shp")
crop(xp, shp)

but I get error that the extents do not overlap! so I try to change the projection of the shapefile as follow, but still I get the same error.

library(rgdal)
library(maptools)
P4S.latlon <- CRS("+proj=utm +zone=16 +datum=WGS84")
indy.shp <- readShapePoly("C:/shapefile/tl_2013_18_cousub.shp", verbose=TRUE, proj4string=P4S.latlon)
writeOGR(indy.shp, dsn = getwd(), layer = "Indy_Bdry", driver="ESRI Shapefile")

I also tried this one, and I get error of "unable to find an inherited method for function ‘spTransform’ for signature ‘"RasterLayer", "CRS""

library(raster)
library(Rgdal)
xp <- readOGR("c:/test/M145.tif")
sh <- readOGR("c:/shapefile/tl_2013_18_cousub.shp")
myExtent <- spTRansform(xp,sh, CRS("+proj=utm +zone=16 +datum=WGS84"))

The link is the shapefile and example of image!

shapefile

image


texts: Openlayers 3 polygon-feature "collision"?

texts:
I am currently building a webapp using openlayers 3.

I build in "GPS-tracking" as a features and I got some questions with that.
So, is it possible to run a function if a GPS-Position as a point-feature, reaches or "collides" with a polygon, which is drawn on the map?

E.g

if(positionFeature.getCoorinates() === polygon1.getCoordinates()){
   return "do someting";
} else {
   return false;
}


texts:
I am learning multithreading and thinking about using it with Arcpy. The example I am working on doesn't have much practical reasons. However, I don't know why it fails on the line "arcpy.Exists(fcPath):" even though the SDE is running and the feature class in question is accessible in ArcCatalog. It didn't help if I used env.workspace.

import threading
import arcpy
from arcpy import env

def tf(conn_string):
    global tTotal
    #env.workspace = conn_string
    fcName = r"GIS_USER.TEST_FC"
    fcPath = conn_string + "\\" + fcName
    if arcpy.Exists(fcPath):
        threadLock.acquire()
        tTotal += 1
        threadLock.release()
    else:
        print('%d: failed to connect to %s'%(tTotal, conn_string))

def mt(conn_string, step):
    global tTotal
    thread_list = []

    for num in range(step):
        t = threading.Thread(target=tf, args=(conn_string,))
        thread_list.append(t)
        t.start()

    for t in thread_list:
        t.join()

    print "total = %d"%tTotal

tTotal = 0
conn_string = r"C:\test_conn.sde"
threadLock = threading.Lock()
mt(conn_string, 5)


texts:
I am using the WMSGetFeatureInfo and I am also obtaining the Info of the base raster and i do not need that. Is it possible to obtain the FeatureInfo of only selected WMS features?

How can I modify this code to accommodate only selected WMS layer?

info = new OpenLayers.Control.WMSGetFeatureInfo({
            url : 'http://niles.iplantcollaborative.org:8080/geoserver/wms',
            title : 'Identify features by clicking',
            queryVisible : true,
            eventListeners : {
            getfeatureinfo : function(event) {
            map.addPopup(new OpenLayers.Popup.FramedCloud("chicken", map.getLonLatFromPixel(event.xy), null, event.text, null, true));
            }
            }
            });
            map.addControl(info);
            info.activate();

UPDATE 2.0

I have added the code,

gfLayers = [];
            var lyr = map.layers;

            for(var a=0; a < lyr.length; a++){
                if(lyr[a].CLASS_NAME == "OpenLayers.Layer.WMS" && lyr[a].visibility == true){
                    if(lyr[a].isBaseLayer == false){
                        gfLayers.push(lyr[a].name);
                    }

                }
            };

            info = {
                hover : new OpenLayers.Control.WMSGetFeatureInfo({
                    url : 'http://niles.iplantcollaborative.org:8080/geoserver/wms',
                    title : 'Identify features by clicking',
                    //layers : gfLayers,
                    layers : [pima, carib],
                    queryVisible : true,
                    hover: true,
                    eventListeners : {
                        getfeatureinfo : function(event) {
                            deletePopUp();
                            popup = new OpenLayers.Popup.FramedCloud("popInfo", map.getLonLatFromPixel(event.xy), null, event.text, null, true);
                            map.addPopup(popup);
                        }
                    }
                })
            };
            map.addControl(info.hover);
            info.hover.activate();

Now I am getting Empty PopUps. not sure why?


texts:
I have 2 shape files,  one that shows all the counties on an island, the other shows all the villages on the island.
Here is the scenario:
My initial map shows the counties. When you click on a county, the map zooms to that county displaying the villages within.

What is the best way to accomplish this using leaflet?


texts: Adding local map as background to dbeaver spatial viewer in offline computer

texts:
I'm using DBeaver (which is a universal database manager) to connect to  a Postgis database on a offline computer (not connected to internet).

The spatial viewer of DBeaver uses Leaflet to display spatial geometries and is very very  useful for debugging purposes. It features the geometries on top of a background OSM map.

My problem: As the computer is not connected to the internet, I have no background map (as there is no connection to OSM).

My Question: Is there a way to config a local background map ? Is there a way I can "make" leaflet go to a local map? (As it is a program, I have no code I can adjust, as far as I can tell).

(I'm using DBeaver version 6.1.5 and Postgis 2.5.2)


texts: ArcGIS Server output directory to mime or not to mime

texts:
Possible Duplicate:
  What are the consequences of setting the SupportedImageReturnTypes to MIME in ArcGIS Server

What are the downsides of changing a map service to MIME only? The developer summit perf video says that it is faster and scales better just using MIME.

What things break from not having the image written to disk if anything?

List of broken things:

WMS
WFS
WCS

I'm also noticing when setting layer definitions the js api posts if the get request is too long and tries to send back the url instead of the mime type. if it's mime only the url json prop is empty.

With arcobjects, if you try doing IMapServerLayout.ExportLayout the url property isn't set but the mimeData property is set.


texts: Creating lines from lonely nodes in a linestring table to closest point on another linestring table

texts:
I want to segmentate a country (denmmark) by using barriers such as larger roads, railroads and the border/ coastline as the delimiting lines of the polygonise function in QGIS.
But as most roads and railroads ends well before hitting the coastline I get some very long polygons along the coast.

I would like to split these long polygons when using the polygonise function by extending road segments close to the border out to the coast. But I only want to select start and end nodes of segments that don't share more than one start and end node with the road network. I've written some SQL for PostGIS that is supposed to solve the problem by utilising that the segments share geometry in the start and end nodes, unless they fulfill my requirements:

create table snaplines3 as
select st_makeline(st_endpoint(a.geom),st_closestpoint(b.geom,st_endpoint(a.geom)))
from vej_jernbane_multi_to_single a,dk_graense b
where st_dwithin(st_endpoint(a.geom),b.geom,50)
and st_touches(a.geom,b.geom) is false
and st_asbinary(st_endpoint(a.geom))!=st_asbinary(st_startpoint(a.geom))
union
select st_makeline(st_startpoint(a.geom),st_closestpoint(b.geom,st_startpoint(a.geom)))
from vej_jernbane_multi_to_single a,dk_graense b
where st_dwithin(st_startpoint(a.geom),b.geom,50)
and st_touches(a.geom,b.geom) is false
and st_asbinary(st_endpoint(a.geom))!=st_asbinary(st_startpoint(a.geom))

Unfurtunately my code doesn't work as the part    st_asbinary(st_endpoint(a.geom))!=st_asbinary(st_startpoint(a.geom)) doesn't seem to do anything.

Does anybody have an idea as to why it doesn't work and alternatively another solution.


texts:
I am using geoserver 2.5.5. I have to get the legends of each category from different link. Now geoserver gives whole legend as a PNG image format using this link http://localhost:8080/geoserver/wms?REQUEST=GetLegendGraphic&VERSION=1.0.0&FORMAT=image/png&WIDTH=20&HEIGHT=20&LAYER=topp:states

I need to have a link which gives each category of legend as a separate image like I need only <2M category in one link like this


texts: How can I get public transport maps in WKT format?

texts:
I can get street maps from openstreetmaps.org (XML) and convert it into WKT with the tool osm2wkt.

There are several people creating specialist public transport maps:

The Transport Map layer on OSM's front page by Andy Allan
Öpnvkarte (öpnvkarte.de) (Also available at openbusmap.org for those with no Ö accessible}
Openptmap (openptmap.org)
OpenStreetBrowser (openstreetbrowser.org)
OSMTransport (3liz.fr) - as per city javascript overlays
LatLon.org's public transport layer
OpenMap.lt public transport layer
OSM Transport Karte Map for Trucks and Cars

However, there is no way to get the XML or WKT format map. How do I get public transport maps in WKT format?


texts: How to merge (mosaic) images that have some similar parts in their file names in R


texts: How to change the encoding once a spatialpolygon is loaded in R and how to speed up loading a geodatabase .gdb file?

texts:
Question
How to change the encoding setting for a spatialpolygon loaded in R without reloading it?  Similar to the argument encoding=UTF-8 when loading the file with readOGR or similar to the drop-down setting in QGIS to change the encoding that QGIS uses to read the file. (I don't need to change the encoding of the file itself as it is saved, just able to read the characters correctly)
Also, what is the most efficient (fast) package/command to load large geodatabase .gdb files?
Details
I have a very large geodatabase file that I want to read in R. Available here: https://www.donneesquebec.ca/recherche/fr/dataset/adresses-quebec, although I don't think the file is a problem. The documentation for the file states that it was in Latin-1 ("Le jeu de caractères sélectionné (le « code page ») est le Windows Latin-1 (windows 1252)"). However, the characters did not load properly with latin1(RÃ©gionale vs. Régionale) and it turned out that UTF-8 was the correct encoding. I had to reload the file each time to try a different encoding and given that the file is large, it was problematic because it look a long time.
#load layer from geodatabase with latin-1
readOGR(dsn=fgbd, layer="Reseau_routier", encoding="latin1", use_iconv=FALSE)
#load layer from geodatabase with
readOGR(dsn=fgbd, layer="Reseau_routier", encoding="UTF-8", use_iconv=FALSE)


texts:
Does anyone have experience with Esri's js charting library called cedar.js?

It is apparently based on D3.

I may be using it for work, but wondering if there are any pitfalls I should know about vs using D3 itself.

Idea is to have dynamic charts as widgets in Portal.


texts: VRT Files from QGIS2+ makes QGIS3 crash when opened

texts:
We got a bunch of VRT files that we made using QGIS 2 and which reference a varying number of ECW files (from 8 to 300).

Those VRT files were perfectly runing QGIS 2.18.

Problem is: Opening those VRT files with QGis3 make either QGIS crash, or they are extremely slow to open (20min+), while QGIS memory consumption grow up to several GB.

We tried to recreate those VRT using the QGIS Build a virtual raster which calls directly gdalbuildvrt, with no success.

Here is the question:

Why does this happen and could it be fixed?
If not, is there a good alternative to VRT?

  NB: An issue already exists at this subject, but no solution was provided (the gdalinfo -approx_stats one does not work).

Here are the infos about one of those VRT:

$ gdalinfo myfile.vrt
Driver: VRT/Virtual Raster
Files: myfile.vrt
       //path/to/Ortho_HR_1342_8200_2015_RVB_CC49.ecw
       //path/to/Ortho_HR_1342_8202_2015_RVB_CC49.ecw
       //path/to/Ortho_HR_1342_8290_2015_RVB_CC49.ecw
      (...plus 300 others...)
Size is 438830, 731382
Coordinate System is `'
Origin = (1341999.944345191121101,8294000.062092189677060)
Pixel Size = (0.199963874261849,-0.199963874081155)
Corner Coordinates:
Upper Left  ( 1341999.944, 8294000.062)
Lower Left  ( 1341999.944, 8147750.084)
Upper Right ( 1429750.091, 8294000.062)
Lower Right ( 1429750.091, 8147750.084)
Center      ( 1385875.018, 8220875.073)
Band 1 Block=128x128 Type=Byte, ColorInterp=Red
Band 2 Block=128x128 Type=Byte, ColorInterp=Green
Band 3 Block=128x128 Type=Byte, ColorInterp=Blue

The content of this VRT:

<VRTDataset rasterXSize="438830" rasterYSize="731382">
  <GeoTransform>  1.3419999443451911e+06,  1.9996387426184947e-01,  0.0000000000000000e+00,  8.2940000620921897e+06,  0.0000000000000000e+00, -1.9996387408115546e-01</GeoTransform>
  <VRTRasterBand dataType="Byte" band="1">
    <ColorInterp>Red</ColorInterp>
    <SimpleSource resampling="nearest">
      <SourceFilename relativeToVRT="0">//path/to/Ortho_HR_1342_8200_2015_RVB_CC49.ecw</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="15003" RasterYSize="25005" DataType="Byte" BlockXSize="256" BlockYSize="256" />
      <SrcRect xOff="0" yOff="0" xSize="15003" ySize="25005" />
      <DstRect xOff="0" yOff="470085" xSize="15003" ySize="25005" />
    </SimpleSource>
    <SimpleSource resampling="nearest">
       (... and so on ...)

  </VRTRasterBand>
</VRTDataset>

  Our configurations are:

  Windows 10 (64bits)
  QGIS 3.4 (64bits)
  QGIS project's csr are RGF93/CC49 (3949)


texts: error:Required argument (index) missing to function: List.get(list, index)

texts:
I have written a code in GEE retrieves images from sentinel 1 and suppose to display those images.
I have created a for loop that suppose to display those images but for some reason I am getting the error message:

  Required argument (index) missing to function: List.get(list, index)

  Returns the element at the specified position in list.  A negative
  index counts backwards from the end of the list.

  Args:

  list (List): Undocumented.   index (Integer): Undocumented.

  Returns the element at the specified position in list.  A negative
  index counts backwards from the end of the list.

I saw that this problem is common but I don't find the problem in my code that led to this error.

this is my code:

//var geometry=geometry2;

// Load the Sentinel-1 ImageCollection.
var sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD');

// Filter by metadata properties.
var vh = sentinel1
  // Filter to get images VH dual polarization.
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
  // Filter to get images collected in interferometric wide swath mode.
  .filter(ee.Filter.eq('instrumentMode', 'IW'))
    // reduce to VH polarization
  .select('VH')
   // filter 10m resolution
  .filter(ee.Filter.eq('resolution_meters', 10))
  .filterBounds(geometry);

  // Filter to get images from different look angles.
var vhDescending = vh.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));
//var vhAscending = vh.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'));

//fiter date
var vhDesc = vhDescending.filterDate(ee.Date('2015-08-01'), ee.Date('2019-10-08'));
//var vhAsce = vhAscending.filterDate(ee.Date('2019-06-01'), ee.Date('2019-10-08'));

//clip the dataset according to the geometry
var clippedCol=vhDesc.map(function(im){
  return im.clip(geometry);
});

//count how many images we have in the collection-
var count=clippedCol.size();
print('Number of Images',count);

//create list of images
var ListDesc =(clippedCol.toList(clippedCol.size()));

var ListofNumbers=[0,1,2,3,5,17,20,28];

//display the images using loop
for (var i in ListofNumbers) {
  var image = ee.Image(ListDesc.get(ListDesc[i]));
  Map.addLayer(image,{}, i);

}

code update:

//var geometry=geometry2;

// Load the Sentinel-1 ImageCollection.
var sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD');

// Filter by metadata properties.
var vh = sentinel1
  // Filter to get images VH dual polarization.
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
  // Filter to get images collected in interferometric wide swath mode.
  .filter(ee.Filter.eq('instrumentMode', 'IW'))
    // reduce to VH polarization
  .select('VH')
   // filter 10m resolution
  .filter(ee.Filter.eq('resolution_meters', 10))
  .filterBounds(geometry);

  // Filter to get images from different look angles.
var vhDescending = vh.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));
//var vhAscending = vh.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'));

//fiter date
var vhDesc = vhDescending.filterDate(ee.Date('2015-08-01'), ee.Date('2019-10-08'));
//var vhAsce = vhAscending.filterDate(ee.Date('2019-06-01'), ee.Date('2019-10-08'));

//clip the dataset according to the geometry
var clippedCol=vhDesc.map(function(im){
  return im.clip(geometry);
});

//count how many images we have in the collection-
var count=clippedCol.size();
print('Number of Images',count);

//create list of images
var ListDesc =(clippedCol.toList(clippedCol.size()));

var ListofNumbers=[0,1,2,3,5,17,20,28];

//display the images using loop
for (var i in ListofNumbers) {
  var image = ee.Image(ListDesc.get(ListofNumbers,i));
  Map.addLayer(image,{}, i);}


texts: Creating a split proportional/graduated symbol map in QGIS

texts:
I am currently trying to make a split graduated/proportional symbol map on Classical music composers by country and musical era and I am having trouble with making the symbol in QGIS. I had to make a .csv table consisting of the totals by country and by era, and I am able to join my this table to a shape file of European countries, but the desired column (totals, representing the total composers of each country from the entire span of all musical eras) from the joined table I want to show represent with the graduated symbol does not show up in the drop down bar on either program.
I am able to get the split symbol using the pie graph feature on QGIS, but I still would like to be able to represent size changes of each pie graph showing the total composers from all eras.


texts: how can i correct atmospherically landsat 8 images?

texts:
I am trying to correct atmospherically an image from Landsat 8 with 6s model from grass 7 and 6.4. Because there are no parameters for Landsat 8 yet, I am trying to give parameters from Landsat 7 images. When I use i.atcorr I get this message:

wavelength  less  than  0.25  micron:
 let's take s(l)=s(0.25)

and the command ends.

How can I correct my Landsat 8 images? Is there someone who knows when Landsat 8 images will officially be in 6s model?


texts:
I need to replace a 'Global Imagery' basemap in Google Maps.

However if I replace it with "Openlayers.Layer.Google", that would use map projection of EPSG:3857 (900913). But I'm working with EPSG: 4326.

This is how the Global Imagery is defined:

       layers: [

new OpenLayers.Layer.WMS("Global Imagery (OpenGeo) [Online]",
                        "http://maps.opengeo.org/geowebcache/service/wms", {
                            layers: "bluemarble"
                        },{

                            isBaseLayer: true,
                            visibility: false
                        }
                    )

]

Is there a way to replace it into Google Maps with EPSG:4326, while keeping the same definition structure?


texts: GeoPandas sjoin largest filter by overlap

texts:
I am using GeoPandas to join two geo dataframes with the parameter how="inner", op="intersects".
Due to nature of the files a larger percentage of the files are duplicates (in regard of the geometry). Therefore I want to ask how to filter out and restrict to a single area, like selecting ownly output file entries with the highest Intersection over Union, DICE coeficient, Overlap coeficient or similar of the overlapping geometries.


texts: Making a tool from script, arcpy.Project_management () comes with an error

texts:
I have a script (just a fragment here, up until the line after which the error comes up):

import arcpy
from arcpy import env

arcpy.env.workspace = arcpy.GetParameterAsText(0)
arcpy.env.overwriteOutput = True

#specify input, output datasets and output coordinate system
input_fc = arcpy.GetParameterAsText(1)
output_fc = arcpy.SetParameterAsText(2, 'ft_mb_afterIteration_lambert')
out_coordinate_system = arcpy.SpatialReference(3308)

#add length column to input datasets
fieldName = "line_length"
arcpy.AddField_management(input_fc, fieldName, "Double")

#project to Lambert and save as new file (for in feature should be Feature Layer / Feature Dataset )
arcpy.Project_management(input_fc, output_fc, out_coordinate_system)

The error is:

Traceback (most recent call last):
  File "H:\Scripts\python\RemoveShortLines.py", line 25, in <module>
    arcpy.Project_management(input_fc, output_fc, out_coordinate_system)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", line 8825, in

  Project
          raise e
      ExecuteError: Failed to execute. Parameters are not valid.
      ERROR 000735: Output Dataset or Feature Class: Value is required
      Failed to execute (Project).

I have also attached the properties of how I set the tool. I am running this script as a tool.

Essentially at that line where error comes up, the the input dataset (input_fc) is converted to a new dataset (output_fc) with different coordinate system.

Can you suggest why my script is not working?


texts: Sum new values to existing attribute values in PyQGIS?

texts:
I'm trying to sum new values from a list in one loop to the existing attribute values (not in a new column, but in the exiting attribute values itselfs).

def radiusDivider(dist):
  try:
    itlayer = QgsMapLayerRegistry.instance().mapLayersByName('itpoint')[0]
    print ('itpoint-layer found')
    prov = itlayer.dataProvider()
    itlayer.startEditing()
    if prov.fieldNameIndex("score") == -1:
        prov.addAttributes( [ QgsField("score", QVariant.Double)])
        print("field added")
    else:
        print ("field already exists")
    itlayer.updateFields()
    it = itlayer.getFeatures()
    count = 0
    for feat in it:
        itlayer.changeAttributeValue(feat.id(), 2, 0)
        #print ("changing attribute")
        count += 1
    itlayer.commitChanges()
    for radius in range (50, dist, 50):
        addToScore(radius)
  except:
    print ("itpoint layer does not exist")

...
a lot of stuff happens in the "addToScore function", but the only important thing is, that I append a calculated “score“ to the "scorelist" for each point of the "itlayer". I'm trying to add in the following code the calculated current score to the existing score (from eventual previous runs) to get after all needed iterations of calculation the final score. :-)
...

def addToScore(bufDist):
    #...
      scorelist.append(calc)
    #...
    itlayer.startEditing()
    itlayer.updateFields()
    it = itlayer.getFeatures()
    count = 0
    idx = itlayer.fieldNameIndex('score')
    for feat in it:
        oldscore = feat.attributes()[idx]
    print ("old:", oldscore)
    print ("current:", scorelist[count])
    newscore = list(np.asarray(oldscore) + (scorelist))
    print ("new:", newscore)
    #prints above show that it works perfectly, so I would guess in the following line is something wrong?
    for ft in it:
        itlayer.changeAttributeValue(ft.id(), 2, (newscore[count]))
        count += 1
    #print ("changing attribute")
    itlayer.commitChanges()
    itlayer.removeSelection()
    #inputlayer.selectByIds(true, 1)
  except IndexError:
    print ('itpoint does not exist')

Edit 1:

I checked and only the values in "scorelist" have type 'int', the other values from "oldscore" and "newscore" have type 'numpy.64int'...? Is this the reason it doesn't work and how do I solve it?