texts: Multiple color styles in one GeoServer layer

texts:
I have a GeoServer layer coming from a shapefile and I am using it with OpenLayers. I want to put multiple colors in the same layer depending on their attributes. How can I put multiple colors on a GeoServer layer depending on the attributes of the layer?

The result I want is something like this:

The only thing I've tried is to create a style in GeoServer and add it to the layer, but I don't know how to add multiple styles depending on the attributes of the layer.

The code of my layer style:

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor version="1.0.0"
 xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
 xmlns:ogc="http://www.opengis.net/ogc"
 xmlns:xlink="http://www.w3.org/1999/xlink"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <!-- a Named Layer is the basic building block of an SLD document -->
  <NamedLayer>
    <Name>poligono_rojo_verde</Name>
    <UserStyle>
    <!-- Styles can have names, titles and abstracts -->
      <Title>poligono_rojo_verde</Title>
      <Abstract>A sample style that draws a polygon</Abstract>
      <!-- FeatureTypeStyles describe how to render different features -->
      <!-- A FeatureTypeStyle for rendering polygons -->
      <FeatureTypeStyle>
        <Rule>
          <Name>rule1</Name>
          <Title>Gray Polygon with Black Outline</Title>
          <Abstract>A polygon with a gray fill and a 1 pixel black outline</Abstract>
          <PolygonSymbolizer>
            <Fill>
              <CssParameter name="fill"> #ff0000 </CssParameter>
            </Fill>
            <Stroke>
              <CssParameter name="stroke">#0cff00</CssParameter>
              <CssParameter name="stroke-width">1</CssParameter>
            </Stroke>
          </PolygonSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts: Why is the length of one degree in Turf 0.9988?

texts:
Calculating the (spherical) distance along one degree of longitude at the equator in Turf returns a number less than 1 (Try it):

require('@turf/turf').distance([0,0], [1, 0], { units: 'degrees' })
0.9988329686371696

Why is this? Is the WGS84 ellipsoid not exactly 360 degrees?


texts:
Writing code using C# and the DotSpatial library, with the following code:

FeatureSet fs = new FeatureSet(FeatureType.Polygon);
        fs = FeatureSet.Open(Application.StartupPath + "\\MapFile\\Land.shp");

I get this error message:

Cannot implicitly convert type 'DotSpatial.Data.IFeatureSet' to 'DotSpatial.Data.FeatureSet'. An explicit conversion exists (are you missing a cast?)


texts: Find the nearest point on polygon (coastline)

texts:
In python I have a shapely polygon (which defines the coastline of the UK), given a point within the UK I would like to find the closest point on the polygon (coastline) to that point. I have been trying to modify the code posted here, but the resulting nearest point I get is incorrect.  What should I change?
from shapely.geometry import Point
import geopandas as gpd
from cartopy.io import shapereader as shp_r
from shapely.ops import nearest_points
import pandas as pd
import geoplot as gplt
import matplotlib.pyplot as plt

# Random point in the U.K (close to the south coast, Worthing)
long, lat = 50.907262, -0.420036
p0 = Point(long, lat)

#10m resolution country data
resolution = '10m'
category = 'cultural'
name = 'admin_0_countries'
shpfilename = shp_r.natural_earth(resolution, category, name)
world_pd = gpd.read_file(shpfilename)

# Filter out United Kingdom
country_name = 'United Kingdom'
country_pd = world_pd[world_pd['ADMIN'] == country_name]

# Extract mainland as single Polygon
uk_main_landmass = country_pd['geometry'].iloc[0].geoms[1]
print (type(uk_main_landmass)) #Out: shapely.geometry.polygon.Polygon

#Plot polygon (UK)
m = gpd.GeoSeries(uk_main_landmass.boundary)
# m.plot()

#Calculate the nearest point on the polygon to the input point (p0)
p1, p0 = nearest_points(uk_main_landmass.boundary, p0)
print (list(p1.coords)) #Result is near Dover rather than on South coast?

nearest_lat, nearest_long   = (p1.coords)[0] #Long and lat are swapped?

#Create geopandas dataframes for random starting point
df_start = pd.DataFrame(
    {'Latitude': [long],
     'Longitude': [lat]})
start = gpd.GeoDataFrame(
    df_start, geometry=gpd.points_from_xy(df_start.Longitude, df_start.Latitude))

#Create geopandas dataframes for nearest point on polygon (coastline)
df_nearest = pd.DataFrame(
    {'Latitude': [nearest_long],
     'Longitude': [nearest_lat]})
nearest = gpd.GeoDataFrame(
    df_nearest, geometry=gpd.points_from_xy(df_nearest.Longitude, df_nearest.Latitude))

#Shows input point (red) and "closest point" (blue) on coastline
ax = gplt.polyplot(country_pd, figsize=(12, 12))
start.plot(ax = ax, color = 'red')
nearest.plot(ax = ax, color = 'blue')
plt.show()


texts: Strange coordinate after use tippecanoe and Mapbox/pbf

texts:
I am trying to make vector tile base on Mapbox pbf standard.

First, convert data from EPSG4326 to EPSG3857 and export as .json.
By using QGIS to export data and set CRS as EPSG3857.

Second, use mapbox/tippecanoe to generate vector tile in pbf format.

Third, use mapbox/pbf, @mapbox/vector-tile and zlib, and use toGeoJSON.

Finally, I got the following results:

geojson {
    type: 'Feature',
    geometry: {
        type: 'LineString',
        coordinates: [[Array], [Array], [Array], [Array], [Array] ]
    },
    properties: {
        FNODE_: 27,
        TNODE_: 27,
        LPOLY_: 9,
        RPOLY_: 1,
        LENGTH: 1.567,
        HYD1_4M_: 278,
        HYD1_4M_ID: 1454,
        GBCODE: 23010,
        NAME: '濂庡悲娌',
        LEVEL_RIVE: 5,
        LEVEL_LAKE: 1
    }
}

and in the coordinates:

[ -48.427734375, 74.74118495548183 ]
[ -48.7353515625, 74.78737860165961 ]
[ -48.6474609375, 74.74118495548183 ]
[ -48.4716796875, 74.71803680830911 ]
[ -48.427734375, 74.74118495548183 ]

This coordinate is very strange, the origin data is in EPSG4326 and I convert the data in EPSG3857, but the above data is obviously not EPSG4326 or EPSG3857. I know the exact location of my data, so this data will not show negative latitude and longitude under EPSG4326.

Does anyone know why this is?


texts: Where deploy Custom Provider DLL for ArcGIS Server if GAC is not an option

texts:
I have a client that has developed a custom Membership Provider for ArcGIS server and cannot deploy the provider to the GAC. Does anyone know where we could put the DLL such that the ArcGIS Security and Token applications can "see" it?

Thanks,
Carl


texts: Any ideas how to stop the Windows version QGIS crashing on exit?

texts:
I’m hoping someone can help.  I know that there’s been a bug with QGIS crashing on exit, but I was hoping that this has been resolved with the latest release of the software (the 64bit version for Windows).  I’ve installed 3.8.0 'Zanzibar' and the LTR 3.4.9 'Madeira', but no such luck with either of these – QGIS continues to crash on exit. I’ve even installed the software using OSGeo4W (64bit), though exactly the same thing happens.

I’ve read somewhere that the best way to come over this issue is to uninstall the version that’s crashing and then to conduct a reinstall.  I thought I had been doing this, but I guess some bogus code is being left on the PC, which could be the culprit...

So, here’s the question: when I’ve done an uninstall do I then need to go into regedit and remove every single occurrence of QGIS and similarly remove the entire contents of the ‘.qgis2’ folder, comprising the cache and python settings etc?

Should any other QGIS folders be removed?

Do I need to completely remove all mention of QGIS from the PC?

Is there anything that can potentially safely stay?

Ultimately, I just want to overcome this issue and I wanted to hear from others on this forum before reporting this.


texts: How do I select * rows of multiple columns in one query using Select by Attributes within ArcMap?

texts:
I am trying to convert a psql query string into arcpy formatting for an UpdateCursor where expression parameter:

'''SELECT \"NAME_0\",\"NAME_1\",\"NAME_2\",\"NAME_3\",\"NAME_4\",\"NAME_5\",\"WIKI_URL\", \"LEVEL_DEPT\" FROM
AdminBoundaries WHERE \"WIKI_URL\" IN(SELECT \"WIKI_URL\" FROM AdminBoundaries GROUP BY \"WIKI_URL\" HAVING
(COUNT (\"WIKI_URL\") > 1)) ORDER BY \"WIKI_URL\";
'''


texts:
I'm trying to develop code that is capable of cutting out a mosaic from the polygons of a shapefile using PyQGIS.

The script is run but does not produce any results.

from qgis.core import *
from qgis.gui import *
import processing
from processing.core.Processing import Processing
import gdal

VECTOR = QgsVectorLayer('C:/Users/phpor/Desktop/','quadricula.shp','ogr')
RASTER = QgsRasterLayer('C:/Users/phpor/Documents/SPU/OS_01/Georreferenciamento_Mosaico_OS_01/Raster/02_Processo/','Aerofoto_50000_Joinville_Mosaico_Recorte.tif')

RASTER.isValid()

row_info = VECTOR.getFeatures()
for row in row_info:
    row_name = row['indice']
    output_file = "C:/Users/phpor/Desktop/teste/" + "cliped_" + str(row['indice']) + ".tif"
    VECTOR.setSubsetString("indice=" + str(row['indice']))

    parameters = {'INPUT': RASTER,'MASK': VECTOR,'NODATA': 255.0,'ALPHA_BAND': False,'CROP_TO_CUTLINE': True,'KEEP_RESOLUTION': True,'OPTIONS': None,'DATA_TYPE': 0,'OUTPUT': output_file}

    clip = processing.run('gdal:cliprasterbymasklayer', parameters)

    VECTOR.setSubsetString('')

print("done")

Now, with some adjusts, I can create a all rasters with the corrects names, but the iteration mask (polygon selected in each iteration) are not respected. In others words, the output of the mask tool are the same for all files (the original draw of the mask).

from qgis.core import *
from qgis.gui import *
import processing
from processing.core.Processing import Processing
import gdal

VECTOR = QgsVectorLayer('C:/Users/phpor/Desktop/quadricula.shp','quadricula.shp','ogr')
RASTER = QgsRasterLayer('C:/Users/phpor/Desktop/aerofoto.tif','aerofoto.tif')
if RASTER.isValid():
    print("ok")
else:
    print("not ok")

Processing.initialize()

row_info = VECTOR.getFeatures()
for row in row_info:
    row_name = row['indice']
    output_file = "C:/Users/phpor/Desktop/teste/" + "cliped_" + str(row['indice']) + ".tif"
    VECTOR.setSubsetString("indice=" + str(row['indice']))

    parameters = {'INPUT': RASTER,
            'MASK': VECTOR,
            'NODATA': 250.0,
            'ALPHA_BAND': False,
            'CROP_TO_CUTLINE': True,
            'KEEP_RESOLUTION': True,
            'OPTIONS': None,
            'DATA_TYPE': 0,
            'OUTPUT': output_file}

    print(row_name)

    clip = processing.run('gdal:cliprasterbymasklayer', parameters)

    VECTOR.setSubsetString('')

 print("done")


texts:
I am currently looking examples of animated maps. Creating animated maps with the help of this post and similar post checking Animaps site also

I am presenting this map in local school so those students should take interest to learn and create various simple interactive maps.

I am in process to create animated map contains Gardens and amusement parks across the city.


texts: How to evaluate performance of GPS/GNSS unit

texts:
I need to integrate some sort of satellite-based positioning for a mobile robot project. I am at the very beginning - I got 2 u-blox GNSS evaluation boards (these NEO-M8Ts) - I chose an integrated antenna and ground plane to get started faster (no external antenna hassle, ground-plane as RTK-ready). I used the u-center software to do some very simple initial tests both measuring stationary and moving positions - I basically only recorded the position from the sensors to a .ubx log file and then played it back.

I am trying to implement this for my mobile robot and I realize that a smartphone-grade GPS (Samsung J5) gives me better preliminary results than an u-blox eval board  - I wonder why, I guess Android may fuse the IMU and have better readings even with worse antenna?

What is the correct way of evaluating GPS performance in terms of position accuracy/precision? I know this will depend on the location and its sky clearance (number of sats used, multipath errors etc.) but is there a rule of thumb? I guess I can only evaluate in stationary position since I will have no correct "ground truth" for mobile testing?

Also, I am interested in RTK performance, but that is probably a candidate for another question.


texts:
I've been upload my own qgis plugin into this:
https://plugins.qgis.org/plugins/Otoklim/

But the plugins has no public version yet..

I didn't get any message from QGIS and i don't know where to ask..

Is there any steps that i should do before uploading to QGIS repository so the plugins can be public version and able to be downloaded via QGIS plugins menu..?


texts:
I notice several EPSG codes in one WKT. At the first glance it almost is non-intuitive to see so many EPSG codes linked to a single EPSG.
Can someone explain each line of the following WKTs for EPSG:5070 and EPSG:4326?
EPSG: 5070

EPSG: 4326


texts: How do you calculate a column value based on the value of another column?

texts:
I am trying to select data based on whether it starts with certain characters and then update a field. For example if [Asbuilt] starts with "WL" then do something. Here are the two methods I have tried and neither work.

dim hyplink as string
If InStr([Asbuilt],"WL") then
   hyplink = "Water Info\Water AsBuilt Scans\City_Waterlines_WL\" & [Asbuilt]
Elseif
   InStr([Asbuilt], "CWL") then
   hyplink = "Water Info\Water AsBuilt Scans\County_Waterlines_CWL\" & [Asbuilt]
Else
hyplink = [hyperlink]
End If

Dim hyplink As String
Dim asbuilt as string
asbuilt = [Asbuilt]
If asbuilt.startswith("WL") then
   hyplink = "Water Info\Water AsBuilt Scans\City_Waterlines_WL\" & [Asbuilt]
Else
   hyplink = [Hyperlink]
End if


texts: Making Tile Service fast in ArcGIS Server for Silverlight Viewer?

texts:
I am publishing a tile service in ArcGIS Server but when I call it in ArcGIS Silverlight Viewer it took a lot of time to open completely, and panning is very slow.

What should I do to make my app smooth and fast?


texts: I use QgsLayerDefinition.loadLayerDefinition(path, project, rootGroup) in my code. But this method changes ids of loaded layers. I need the codes to be old. Is there any way to change ids (for example any custom) of layers after loading layer definitions? Or have any other ideas?


texts: Is there public OGC Catalog Service (CSW) available which lets ArcGIS for Desktop with CSW Client Add (WMS) To Map?

texts:
This is a follow on from a Question that I asked last year entitled Is there a public OGC Catalog Service (CSW) available?.  That Question has an Accept-ed Answer but I would now like to expand upon the question slightly.

Is there a public OGC Catalog Service (CSW) available which serves URLs that can be bound to as Live Data and Maps by ArcGIS for Desktop 10.2?

I would like to learn of a CSW service URL (with its Profile specified) that I could use to Configure the CSW Client of ArcGIS for Desktop 10.2 so that I can see it add a WMS service to ArcMap by using the Add to Map button rather than by using copy/paste of URLs.

A Comment on my Answer to Is there an exhaustive, searchable catalog of all GIS web services (e.g. OGC WMS/WFS, REST, SOAP)? says that it is possible but I do not recall ever being able to do it personally, including during tests I performed yesterday using the latest versions of ArcGIS for Desktop and its CSW Client.

After upgrading to ArcGIS 10.3 for Desktop and @Luke resolving Where to download CSW Client for ArcMap (ArcGIS for Desktop)? at that new version, I was able to conclude my test successfully.

I cannot comment on whether it now works at 10.2.x because I forgot to re-check before upgrading from 10.2.2 and my initial testing was at 10.2.


texts:
I am trying to remove a join in arcpy and have been all over the forums and help topics but cannot get my script to remove the join.

I have tried many combinations of inputs to the remove join but to no avail.

Name94 = "C:\\a\\b\\c\\input1.shp"
Name95 = "C:\\a\\b\\c\\input2.shp"
TablesView = "C:\\a\\b\\c\\TablesView.lyr"
arcpy.MakeFeatureLayer_management(Name94, TablesView)
arcpy.AddJoin_management(TablesView, "Buf_Units", Name95, "Units", "KEEP_ALL")
arcpy.CalculateField_management(TablesView, "Buf_Area", "[Area]", "VB", "")
arcpy.RemoveJoin_management(TablesView, "Buf_Units")

This yields the error:

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000800: The value is not a member of  | Input2.
Failed to execute (RemoveJoin).

In place of buf_units I have tried:

arcpy.RemoveJoin_management(TablesView, "'"+Name95.rstrip(".shp")+"'")
arcpy.RemoveJoin_management(TablesView)
arcpy.RemoveJoin_management(TablesView, "input2")
arcpy.RemoveJoin_management(TablesView, Name95)
arcpy.RemoveJoin_management(TablesView, "units")

The previous error is the result of most of these tries but when it is not it is:

ExecuteError: ERROR 000229: Cannot open C:\a\b\c\TablesView.lyr
Failed to execute (RemoveJoin).

I am running python 2.7, 64bit, 10.1 arcgis


texts: Changing layer properties for 3D View using PyQGIS

texts:
I would like to change the default settings in 3D View.

No symbols to Single Symbol
Altitude clamping from Relative to Absolute
Change colors for Diffuse, Ambient and Specular

Changing by hand works fine but is impractical for 10+ layers. I would like to automate it with Python.
I have looked into the 3D Library Python API but I have not been able to get this to work.

I have written some code:

import qgis._3d as d

wkt = "MULTIPOLYGON Z(((0 0 0, 0 0 1, 0 1 0, 0 0 0)), ((0 0 0, 0 1 0, 1 1 0, 0 0 0)))"

symbol = d.QgsPolygon3DSymbol()
symbol.setAddBackFaces(False)
symbol.setAltitudeBinding(1)
symbol.setAltitudeClamping(0)
symbol.setCullingMode(0)

renderer = d.QgsVectorLayer3DRenderer()
renderer.setSymbol(symbol)

materialSettings = d.QgsPhongMaterialSettings()
materialSettings.setAmbient(QColor(255,0,0))
materialSettings.setDiffuse(QColor(255,0,0))
materialSettings.setSpecular(QColor(255,0,0))
symbol.setMaterial(materialSettings)

layer = QgsVectorLayer("multipolygon?crs=EPSG:31469", "MULTIPOLYGON Z", "memory")
renderer.setLayer(layer)
feature = QgsFeature()
geometry = QgsGeometry.fromWkt(wkt)
feature.setGeometry(geometry)
layer.dataProvider().addFeatures([feature])
QgsProject.instance().addMapLayers([layer])

It shows no errors, it changes nothing.

I am using QGIS 3.12.0 Python console.


texts: Controlling Image Output Settings in QGIS?

texts:
I am exporting a QGIS map to an image format. I would like to have control over the output settings of the image (dpi), as currently it is only exporting the images at 96 dpi.

Does anyone know if this is possible?


texts: GDAL: using the -R function in raster2psql creates errors

texts:
I'm loading jp2 rasters into my PostgreSQL/PostGIS, and everything works fine unless I use -R in the commands:

raster2pgsql -s 28356 -I -C -M -R -l 8 10cmimagery.jp2 -F -t 250x250 imagery.10cmimageryl8r | psql -h localhost -U postgres -d GIS

The error is:

ERROR:  rt_band_load_offline_data: Cannot open offline raster: 10cmimagery.jp2
CONTEXT:  SQL function "st_transform" statement 1

If I omit the -R there are no errors and the jp2 file imports correctly. I'd like to read the raster from the original file rather than import it into the database ideally.


texts: Setting minZoom at runtime in ArcGIS API for JavaScript

texts:
I want to set maxZoom and minZoom at runtime.

Initially I set maxZoom and minZoom as shown in below code

new this.Map("map1", {
    basemap: "streets",
    center: new this.Point3(-98.35, 39.50),
    zoom: 5,
    maxZoom: 20,
    minZoom: 5
});

After some time user want to change the Min and Max zoom levels at runtime.
So I tried as below code, but not workout.

this.selMapObj.setView(new this.MapView({
    maxZoom: 20,
    minZoom: 9,
}));


texts: Public GIS data for lots within sections of township and range for NM, USA

texts:
Anyone know where to find a publicly available shapefile/fc of lot divisions within sections for township and range of New Mexico?

I've been to several government websites that offer this type of data but can't seem to find any files that are broken down all of the way to the lot level.  I.E. land grids for section and twp/rng seem easy enough to come by from these sites where as lots are eluding me.


texts: Very poor performance of QGIS and PostGIS 2. User error or feature?

texts:
I have a large table in PostGIS 2. It has >500 million polygon, line and point features all happily living in there and the table performs brilliantly in other GIS.

I am now trying to just add that table to a project in QGIS but here is the thing, why does it take so long?

When I try to 'Add PostGIS overlay' it says detecting for about four hours to tell there are three geometry types.  OK, I can take that if I only have to do it once...

Once it finishes that I select all three parts click OK...  So far I'm two days into trying to get it to load into the canvas.

What is going on?  I am zoomed in so only a few polygons should be rendered on screen.

Is this an error or a 'feature' of QGIS?  Is there something I need to specifically do to enhance the performance for QGIS?

I am using MapInfo and Cadcorp SIS both of which love the "everything in one place" model rather than the splitting it up, join it back later with queries model.  And thanks @AndrewJoost, the estimated metadata saved about an hour but still it takes a few hours too long and misses one of the geometry types off the list...   QGIS not supporting mixed geometries in one overlay is no problem, as long as it supports them in one table.

And the way it does this with a filter is great for small tables but just fails when you have a large one it seems.

I am unable to separate the table into three parts as this does not make sense for the other systems using the data, and would also break them.

@NathanW How can I get 3 records for a single table into the geometry_columns view now that it is no longer a physical table? any ideas?  Sorry about all the newbie questions

Any other ideas to stop QGIS trying to be helpful and scanning a table before allowing me to add it? I don't mind choosing geometry types by hand for every overlay if that's what it takes?


texts: QGIS Layout Print Composer: getting the height (etc.) variable of a certain layout item

texts:
Currently I am drawing two polygons that are supposed to be the "shadow" of legend since apparently there is no shadow on frame option.
Yes this question is related to this one QGIS print layout: creating a shadow arroud a frame for legend at only specific sides, but considering the topic is different I am asking a new one here.
For now I am just testing it with one.
The problem is however that the height and width of the Legend is dynamic depending on what gets into it, so to sure the polygons also adapt to the legend, I will be giving them the height and width values of the legend plus a certain marge. Now in theory I should be using this formula because my Item Id of the Legend is Legend.
map_get(item_variables( 'Legend' ), 'height')

However the problem is that this does not actually work. Not to mention I need to add some extra values like 0.3 (because of the frame size of the legend) so in theory I should also do this:
map_get(item_variables( 'Legend' ), 'height') + 0.3

What am I doing wrong here because when the legend auto adapts the polygon does not adapt?


texts: Calculate the area covered by a pixel in Google Earth Engine?

texts:
I am doing a project on calculating the growth of urban areas in different cities located in different parts of India.

I have been able to classify the area into urban, vegetation and water.

I cannot calculate the area covered by each pixel. I searched on Google and the available methods use band values, what I need is to find the area covered by urban/water/vegetation based on spectral end-member values that are been defined.

If anyone has the idea on how to perform it please share and if possible suggest an idea to export it into a graph.

Here is a link to the script:

https://code.earthengine.google.com/d5e7c1435a160524404878c6e2110bda

EDIT: my aim to find the growth of an urban area or loss of vegetation in a region over the years. The urban area is indicated in red colour when you run the script.

I am finding how much area is covered by urban area so that I could measure the growth of the region.

I'm fairly new to Earth Engine.


texts:
I'm trying to convert one raster into a SpatialPolygons. I've got a big raster (12 millions pixels) with a lot of zeros, and some negative values (representing some depression in a DEM). I would like to cluster all these negatives values into polygons, to get a map of all the depressions on my map.

I'm a R beginner and I can't find a suitable function to do the job.

Does anyone have a clue ?


texts: Leaflet printing does not show pins and map is shifted

texts:
I don't exactly why but whenever I try to print my leaflet map with layers, it shows the map a little shifted to the left as well as the pins are not quite ok. I tried using easyPrint() and window.print(), same result. It used to work well. Any ideas? Here is the difference

when it opens up to print the map


texts: How can I calculate the distances between my points in sequence automatically?

texts:
I am trying to measure distances between points and writing the calculated measure between these points in the attribute table.

The points represents a vehicle's location based on GPS data according to existence location in time aspect. I could add the longitude and latitude data from Excel to a shape layer.

After creating shapefile, I could be able to measure the euclidean distances between the nodes, but if I will do the calculations manually, I might need to have immortal life.

How can I calculate automatically the distances between the nodes and write the calculations in my table or Excel or CSV data.

the distance will calculate according to the id's.

I mean the first distance will be for FID 1 and 2; the second will be for 2-3, third 3-4 and etc.

Please help me :)


texts: Extract an "average polygon" from multiple right-turn loop polygons?

texts:
Now we are using OSM and OSMnx python package to do a job of extracting the right turn loop in a designated area. We have extracted all the right-turn loops features in the area(Such as the number of nodes, angle, edge length).
We want to extract a polygon shape that represents the characteristics of the whole right turn loops in this area. But these right-turn loops have different polygonal shapes, such as triangles, quadrilaterals, etc., it is not possible to directly handle all right-turn loops.
One idea is to classify these right-turning loop polygons according to the number of nodes before processing. But I encountered the following problems:

At first I want to directly average the included angle and edge length of polygons with the same number of nodes to get the final polygon. But there are many curved roads in the actual roads, which makes the included angle meaningless.

How to get a representative right-turn loop polygon from these right-turn loops using these features(nodes number, included angle, edge length)?

Maybe we can’t solve the problem by using these features. Is there any other better way to solve this problem?


texts: Clipping raster by value of overlapping raster in QGIS?

texts:
I have two overlapping raster dataset with the same resolution. First one is classified into 1 and 0 values. I would like to clip the second raster based on the pixels of the fist raster with value equal to 1.

Is there any tool that can do this job in QGIS?


texts: What is the convergence threshold for ArcGIS Iso Cluster?

texts: The ArcGIS tool Iso Cluster Unsupervised Classification does not seem to have any documentation on the convergence threshold used in the tool.  To compare, Erdas Imagine has all of the necessary parameters displayed.  How does ArcGIS determine the convergence threshold for their version of the ISODATA algorithm?  More importantly, is there a way to set or override the convergence threshold in ArcGIS?


texts: Why are my area calculation values so low?

texts:
My CRS is US national Atlas Equal Area | EPSG:2163

First, I right-clicked on my layer, opened the Attribute Table, toggled editing mode, then opened the Field Calculator.

Next, I checked the box for Create a new field, put in "area_calc" for the Output field name, then chose $area as the value.

However, the problem I'm having is that the area calculation values don't really make sense (I think). For example, the value for area of Alabama on Wikipedia is 135,765 km^2, and the area_calc value is 13. All fields in this column are extremely low compared to the actual state area values.

My question is:

1) If the CRS units is in meters, why are the area values so low?

NOTE: When I use Vector -> Geometry -> Export/add Geometry Columns, the results are the same.


texts: How to rasterize individual feature (polygon) from shapefile using GDAL/OGR in Python?

texts:
I would like to convert the individual records (polygons/features) of a shapefile to raster using GDAL/OGR. Here is the example shape and the mask raster. Rasterizing the full layer with all records (two, in this example) at once works:

from osgeo import gdal
from osgeo import ogr
from osgeo import gdalconst

raster_file = r"c:\Temp\mask.tif"
shape_file = r"c:\Temp\test2.shp"

# open raster and obtain extent and properties
data = gdal.Open(raster_file, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * data.RasterXSize
y_min = y_max + geo_transform[5] * data.RasterYSize
x_res = data.RasterXSize
y_res = data.RasterYSize
pixel_width = geo_transform[1]

# open shapefile and create layer
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shape_file, 0)
layer = dataSource.GetLayer()

# convert the layer with field 'presence' to raster
target_ds = gdal.GetDriverByName('GTiff').Create(r"C:\Temp\test2_rasterized.tif", x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))
gdal.RasterizeLayer(target_ds, [1], layer, options=["ATTRIBUTE=presence"])

When I loop over the individual features, I can't find a way to export the feature to raster:

from osgeo import gdal
from osgeo import ogr
from osgeo import gdalconst
import os

raster_file = r"c:\Temp\mask.tif"
shape_file = r"c:\Temp\test2.shp"

# open raster and obtain extent and properties
data = gdal.Open(raster_file, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * data.RasterXSize
y_min = y_max + geo_transform[5] * data.RasterYSize
x_res = data.RasterXSize
y_res = data.RasterYSize
pixel_width = geo_transform[1]

# open shapefile and create layer
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shape_file, 0)
layer = dataSource.GetLayer()
nr_features = layer.GetFeatureCount()

for feat_id in range(nr_features):
    feature = layer.GetFeature(feat_id)
    output = os.path.join(r"C:\Temp", "rasterized"+str(feat_id)+".tif")
    target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Byte)
    target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))

    # this does not work since I do not pass a layer to the command
    gdal.RasterizeLayer(target_ds, [1], feature, options=["ATTRIBUTE=presence"])

    # e.g. is it possible to achieve it with Rasterize?
    gdal.Rasterize(...)

Is what I am after at all possible using GDAL/OGR?
Or do I need a workaround like converting each feature to an individual layer and then rasterize each layer?


texts:
In the QGIS version 2.X, there was an Iterative Button from the algorithm also available for the algorithm used inside the processing modeller.

Actually, I use version 3.2 and 3.4 and can not find it. Will I find the button somewhere else or is the iterative button removed?


texts: I have an OpenStreetMap PostGIS database and would like to select specific data (for example all shops) of each tile (on a defined zoomlevel) on the planet. Is this possible to achieve and if yes, how?


texts: How to convert Shapefile to KML in Python

texts:
My goal is to extract certain shapes from a shapefile, and render the chosen shapes in a KML file. From searching gis.stackexchange, I was able to piece together the necessary components, and the answer I came to seems to be very close. However, a shape in Gilbert, Arizona is rendered in Kogi, Nigeria. Are the ArcGIS point coordinates relative to a specific reference point?

For the shape files themselves, you can see my project at https://github.com/pmains/geo_shapes

    from pyproj import Proj, transform
    import shapefile
    import simplekml

    def esri_to_geo(*points):
        # ArcGIS format
        in_proj = Proj(init='epsg:3857')
        # Standard geodetic format
        out_proj = Proj(init='epsg:4326')
        # New array that will hold points in latitude/longitude
        new_points = []
        for point in points:
            x2, y2 = transform(in_proj, out_proj, point[0], point[1])
            # Append points as tuples
            new_points.append((x2, y2))
        return new_points

    def main():
        # Read in all precincts
        sf = shapefile.Reader("VotingPrecincts.shp")
        shape_records = sf.shapeRecords()

        # Find Gilbert Gardens (Precinct 251)
        for i in xrange(len(shape_records)):
            sr = shape_records[i]
            if 'GILBERT GARDENS' in sr.record:
                break

        # Create KML file
        kml = simplekml.Kml()
        kml.newpolygon()
        poly = kml.features[0]
        # Convert points
        pts = esri_to_geo(*sr.shape.points)
        poly.outerboundaryis = pts
        kml.save("GilbertGardens.kml")

    if __name__ == "__main__":
        main()


texts: Create new network dataset in Arcgis Pro

texts:
I am trying to create a new network dataset in ArcGIS Pro. I have created a feature dataset and imported the roads layer in it. After that i created topology and fixed the issues in road dataset. After that i am trying to create new road network by right clicking the feature dataset and there is no option to do this.
i looked for tools in toolbox the only option i get is to create a network dataset from some template

Where i can get that template?

Why cannot i use my own schema?


texts: CartoDB with KoboToolbox, Open Data Kit (ODK), Enketo, or Ona source data?

texts:
I'm trying to use survey data collected by Kobotoolbox with CartoDB. I've reviewed the API's docs (https://kc.kobotoolbox.org/api/v1/), and I'm doing my best to follow them to connect my survey's data, but I can't seem to get it right.  I'm left wondering if it's even possible.

You can access your data via basic and token authentication.  I would like to do token authentication.  To access your data via token authentication, you have to register it.  .  I've registered my application with the following details:

redirect_uri:  [username].cartodb.com
clientID:  xxxx (40 characters in total actually)
response_type should be "code"
state:  "a random state string that you client application will get
when redirection happens"

The example given for a successful token authentication to a user's protected data in the KoboToolbox API is as follows: http://localhost:8000/o/authorize?client_id=xxxx8&response_type=code&state=xyz.

Also, a successful call to the data for a specific survey looks as follows:
curl -X GET https://kc.kobotoolbox.org/api/v1/data/yyyy (where yyyy is the unique ID for that survey)

If that's the case, should I just be able to insert my call to my specific survey into the "connect" to dataset field in cartoDB in this format?

https://kc.kobotoolbox.org/api/v1/data/yyyy.csv/o/authorize?client_id=xxxx8&response_type=code&state=xyz

If that is not the correct format, what format should i use?  Or am i trying this completely incorrectly, and what should i be doing then to access data from KoboToolbox?

Fyi, KoboToolbox is one of a number of resources in the OpenRosa ecosystem of open source tools for data collection in areas without data collection - this ecosystem includes Open Data Kit (ODK), FormHub (depricated) Enketo, Ona, SMAP server, and others.  If users have experience with one, I imagine the APIs for others would be very similar, and I'd be interested in hearing from anyone with experience working with them and CartoDB.


texts: Using PostGIS 2.1, getting "Could not load the output GDAL driver"

texts:
I have problem of loading GDAL in PostGIS

SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';

After running

ALTER DATABASE "OME"
  SET postgis.gdal_enabled_drivers
  TO 'GTiff PNG JPEG';`

I still have an error:

rt_raster_to_gdal: Could not load the output GDAL driver

Note that command gdalinfo works. Here is the result executing postgis_full_version() in pgadmin:

  POSTGIS="2.1.7 r13414" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0,
  6 March 2012" GDAL="GDAL 1.11.1, released 2014/09/24 GDAL_DATA not
  found" LIBXML="2.7.8" LIBJSON="UNKNOWN" TOPOLOGY RASTER


texts: What is the alternative of setZoomLevel in vector based map?

texts:
It is commonly used to generate map tiles at a certain scale and re-organize them for web map on PC. For example Google map, and Bing map.

Then user will use the setZoomLevel(int level) or something else to set the zoom of the map. This is because the tiles is not continuous.

Now I wonder what is the alternative method when use display the vector data ?


texts: Exporting Field Properties summary report for shapefile or geodatabase layer?

texts:
I have a File Geodatabase from an external source that has a number of layers.  I need to crosswalk this data to our own data model so some or part of it can be imported into our own database.

To do this, I want to first create a "field map" in Excel that has a separate tab for each of their layer attribute tables that lists their field name, the field data type (long, text, double, etc.), and the field length.  I will then add a column and use that to populate which of our own fields are most equivalent for each one (if any).

This would be much faster if there was a way to export just the field names and their properties for each layer from ArcMap as a table - does anyone know if that's possible?

Right now I'm just exporting one data record for each table, opening that in Excel, then copying the field names row and pasting it as a "transpose" to turn it into a vertical list.  I then would have to look up each field one-by-one in the layer properties window and manually enter in the field type and length for each one.  There will probably be well over 100 fields though, so this is going to be a pain doing it that way.


texts: Is there a way to calculate closest distances to a point without using Near Analysis?

texts:
I frequently get asked to calculate distances from a specific point (A drilling well in the middle of the sea for example) to a number of Environmental layers (Protected Areas, Pipelines, Other wells etc, so would include points, lines and polygons). I currently do it manually which takes time so I thought there could be a way to automate this using python. It seems the best way to do this is use Near Analysis, which I do not have licence for.

I was just wondering is there anyway to use python to still generate some sort of table that would give me the closest distances from my input feature to my numerous environmental features, the names of said closest features, and the distance (and also direction if possible, but thats not essential).

Any ideas?


texts: Valid ranges for bounding box in WMS 1.1.1

texts:
I'm connecting to a WMS server and found that the longitudes LatLonBoundingBox specified for some of the values are slightly outside of the +/-180 range (such as 180.05). Are these values considered valid (I couldn't find anything in the schema that seemed to enforcer a range)?

If they are valid, should they be treated as rounding errors (maybe the values were auto generated) and clamped to +/-180, or should they be normalized in the +/-180 range, so 180.05 might become -179.95?

Also note that in this particular case, the parent layer's bounds are listed as -90/90 and -180/180.

Thanks.


texts: Google Spreadsheet as Configuration for Leaflet App

texts:
I would like to build a Leaflet Webapp that is highly adaptable and easy to configure.
Therefore I would like to use Google Spreadsheets to store the configuration and retrieve the stored information during initialization of the Leaflet-Map to set the used layers, initial centerpoint and so on.

Has anyone done this before or could give me hints?
I have googled a lot about it but haven't found anything like this so far.


texts: How to add to an existing dataset in Model Builder?

texts: I want to run a custom tool created with the Model Builder in ArcGIS.  It needs to output to the SAME file each time it runs, just adding the data.  Basically I have a custom Dissolve tool that will dissolve (in this case) counties together based on user selection.  I need a single file containing ALL the dissolved final shapes, but the Model Builder complains unless each output file is unique. Is there anyway around this or do I have to create dissolve_1, dissolve_2, etc... (to 500 or so) and then combine them afterwards?


texts:
I have attempted to build a virtual raster in QGIS to avoid working with several smaller mosaics.

That ended up well, and I saved all of these individual files as tiff files.
Now I am trying to save all of these as one large file, but whenever I attempt to VRT, the result only shows up one of the several images.

How is that?


texts: Speeding up large GeoJSONs in Leaflet maps

texts:
I am trying to use Leaflet to make a map which displays multiple projection boundaries, i.e. UTM, State Plane, MGRS. Right now I have three geojson files inside the javascript, but they are quite large and I want to add even more. I need to find a way to add all of these different layers while maintaining performance and speed. I'm not sure if that involves creating tiles or asynchronously loading the json files; either way I'm not sure how to do that. I could use some advice and/or links to resources and tutorials that may help.

Here is a link to my site and what already exists: http://www.paulcflynn.com/projection_picker.html


texts: How to randomly select from existing points within each polygon?

texts:
I have a polygon layer (where each feature has a class attribute), and a point layer, and there are at least 10 points within each feature of the polygon. Is there a way to randomly select only "k" points of these 10 points for each feature (class) in QGIS?

Thanks!


texts: Why is GDAL ECW compression much poorer than ECW compression from ERDAS

texts:
We use ECW compressor shipped with ERDAS Imagine to export Orthoimagery for our different departments. Because we are facing problems with the Imagine license server, we had to choose GDAL (built with ECW SDK) as a fallback solution. Because the Orthotiles are all approximately 100MB large, thats possible.

As a compression ratio we use 1:20. With ERDAS ECW compressor (it does not matter if we use the licensed one or the free edition from some years ago) the images are all compressed to approximately 4 to 6 MB each.

With GDAL and the same compression ratio, the images got 10 to 20 (!) MB each. ALso some of the files cannot be opened afterwards. Also when I increase the compression ratio (to 1:30 for example) the imagefiles get smaller but face a lack of imagequality very quickly.

Why is GDALs compression significantly more worse than ECW Compressor, when they are both based on the same libraries from the ECW SDK?

Is it because GDALs implementation is so buggy or because ERDAS is keeping its Aces on their sleeve?

EDIT: Yeah, I know that this is truly a naive question ;)


texts: Why does Clip keep giving Error 999999: Error executing Function Invalid Topology [Topoengine Error]?

texts:
The layer being used as the Clip Features has been used before and worked fine, so I assume the issue is the layer I am trying to Clip (Input Features) - I keep receiving the same error when running the Clip tool :

  Error 999999: Error executing Function Invalid Topology [Topoengine
  Error]

The Clip tool has worked fine prior to this. After reading online similar queries regarding this error, the solution in the majority of cases was to change the Input Features to a shapefile - I have done this but receive the same error.

What is my best workaround for this issue?


texts: Saving a vector layer as .xyz with PyQGIS 3

texts:
I'm writing a standalone script for QGIS. In this script I want to import a couple of text files saved as .xyz. When I import them in QGIS I have to use "Import layers as delimited text", thus they are not rasters. They are text files with lots of rows and each row contains x,y and z coordinates.

After I import them I want to check a condition and export them to a new folder.
The import part works fine but I get no export files in my output folder. That's the export part of the script:

lyr = QgsVectorLayer(uri, name, "delimitedtext")
outDir = 'C:/Users/user/Desktop/output'
fields = lyr.fields()
QgsVectorFileWriter.writeAsVectorFormat(lyr, outDir, "utf-8", fields, "CSV",
                                        layerOptions=['GEOMETRY=AS_XYZ'])

When I try to run the script, I don't get an error message but there are also no files in my output folder. I don't see the error. Can someone help me?

I'm using QGIS 3.4.13

EDIT

I checked the answers and modified the code but it is still not working. Still no error message and still no files in my output folder.

This is the whole function without the condition:

def moveFiles(self):
    inputDir = self.dlg.input_folder.text() #directory of xyz files - "Users/denni/Desktop/xyz_data"
    outDir = self.dlg.output_folder.text()  #output directory - "C:/Users/denni/Desktop/output/"
    crs = "EPSG:25832"

    for file in glob.glob("*.xyz"):
        uri = "file:///" + inputDir +"/"+ file + "?type=csv&delimiter=%s&crs=%s&xField=%s&yField=%s" % (" ", crs, "field_1", "field_2")
        name = file.replace('.xyz', '')
        lyr = QgsVectorLayer(uri, name, "delimitedtext")
        outputPath = outDir + file
        QgsVectorFileWriter.writeAsVectorFormat(lyr, outputPath, 'utf-8', lyr.crs(), "CSV",
                                                layerOptions='GEOMETRY=AS_XYZ')

What could be the problem?


texts: Use a field from atlas layer in aggregate function

texts:
I cant find a solution to use the value of the current atlas entity in the aggregate function in the filter expression part.

Example:

aggregate('layer_to_aggregate','count',"field_form_layer_to_aggregate_to_aggregate","field_from_layer_to_aggregate_to_filter"="field_from_atlas_layer")


texts: Can I edit coordinates within an edit session and reproject without reloading the data

texts: I have one lat/long to correct.  I edited within my table, saved edits, and now when I project using the Project tool the point stays on my map in the same position.  Do I need to reload my data from an edited Excel file or can I do it from within Arc?


texts: GeoDjango - return SRID of geometry in GeoQuerySet

texts:
I'm wondering if there is a function in GeoDjango which would allow me to determine the SRID of a geometry column in a GeoQuerySet.
I know I could do it with a raw SQL query (I'm using PostGIS) but I have not been able to find a GeoDjango function which would allow me to get the SRID easily.

Thanks


texts: Interpolating orthodrome between two (lon,lat) points in Python?

texts:
I have a polygon that consists of a few (lon, lat) points, each pair of consequent points is considered to be connected by orthodromes. I want to interpolate the polygon with some extra points between each source point pair (add some more waypoints lying on the orthodrome).

I know, that geopy, for example, has the possibility to calculate the distance on great circle, but I need some more: to find the waypoints on great circle.

I need something like that:

src_lon = 10
src_lat = 40

dst_lon = 20
dst_lat = 60

n_points = 100

lons, lats = interpolate_orthodrome(src_lon, src_lat, dst_lon, dst_lat, n_points)


texts: Accessing geojson properties using OpenLayers3

texts:
Is it possible to access geojson properties using the OpenLayers library? I tried to search for it but couldn't find any proper answer.
According to my research, the library converts each geojson feature to ol.Feature which doesn't provide any access to the properties.

For eg:

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Point",
        "coordinates": [102.0, 0.6]
      },
      "properties": {
        "prop0": "value0"
      }
    }
}

How can I access prop0 using the openlayers library?


texts:
I have a leaflet marker created similar to the example below.

var marker = L.marker([map.getCenter().lat,map.getCenter().lng], {
      alt: tempID,
      icon: starIcon,
      draggable: true
  })
  .addTo(map)

I would like to toggle the draggable property. I've tried accessing manually like below to no avail.. Any suggestions?

marker.options.draggable=false;


texts: What are my options to run a WFS with Geoserver and IIS?

texts:
I want to set up a WFS using Geoserver within Jetty servlet container. The server is running Windows Server 2008 R2 and is using IIS primarily. Not opening another port is preferable. Apparently one solution could be utilizing IIS’s  Url-Rewriting to make Geoserver (or the WFS) accessible. The WFS should only be available for authenticated users – I have to use digest-encryption.

When WFS is working there are basically two uses cases:

FME-Script should use the WFS to retrieve data to process
web-app monitoring the data (Openlayers 3)

In the first case the access should not be the problem. I am a bit worried regarding the authentication. Basically FME can handle digest-encryption, but can the url-rewriting possibly cause problems? Can anyone share experiences?
In the second case I do not know if it is even possible. The web-application is going to be on a completely different server, so the url rewriting cannot serve as a (reverse) proxy. I might use a proxy here or use cross domain access. Is there any other option? How can I handle the authentication – is it even possible?

Another thing that comes into my mind is: Do cookies/session-object still work, or must Geoserver run through the whole authentication chain again for every request?


texts: FME 2015 - how to write attribute name in a field?

texts: I am merging 5 csv files which each contain one field and need to sort them according to the length of the field title. This means I need to first read the attribute name from the csv file and write it to a second file in order to create another field containing the length of that field name. I can´t work out how to expose the attribute name though. Does anyone know how to do this?


texts: Difference in area size calculation since update packages

texts:
On my mac I use brew to install PostGIS and all dependencies, including geos.
Since month or so, I have a failing test in my Python code.
A polygon that is entered in WGS84 (srid 4326), the converted to srid 28992 (RDS new) to get the surface in km2 used to give a surface of 76.33km2. Since updating packages, the surface is 225.14km2. In PostGIS, the calculation is just fine.
SELECT ST_Area(ST_Transform(ST_GeomFromEWKT('SRID=4326;POLYGON ((4 52, 4.1 52, 4.1 52.1, 4 52.1, 4 52))'), 28992))/1000000
# result: 76.32601969674697

In Django, this calculation is "outsourced" to the geos package, and gives 225.14km2.
AssertionError: False is not true : The field 'geom' on form 'form' in context 0 does not contain the error 'Geometry too big: area is 76.33km², maximum area is 25km².' (actual errors: ['Geometry too big: area is 225.14km², maximum area is 25km².'])

I have rebuilt my environments, reinstalled everything with brew, but I'm completely stumped as to how this can happen. in our CI process, the same tests passes (GitHub actions).
I'm looking for anything I can check to see if the projection used by geos on my mac is correct, or otherwise can help me find the reason for this change.


texts: Choosing a Projected CRS in Order of Making a 3D Polygon Display in QGIS 3.4.1

texts:
I have been followed the Instructions explained how to make 3D polygons.

I have downloaded some random polygons from OSM and defined the 'Extrusion' with the expression rand (5, 40).

Then, I select the "3D View" function.

At this point I got an error message that telling me that my CRS are not "projected CRS".

I used in this project "WGS 84 / Pseudo-Mercator" CRS, which I do not know to which CRS I should change them into. I changed the measurement units to degrees, it did not help.

What is the correct CRS?
Is there are anything more need to be done in order of making a 3D map from a random polygon that is not listed in the steps I have listed?


texts: Creating catalog of vector layers using ArcGIS Platform?

texts:
Is there a way to create a catalog of vector layers in ESRI? We have thousands of tiles of 25cm contour data which we need to individually find and add to maps. It would be so much easier if there was a way to load the best available tiles based on the current map extent.

For Raster we use the Mosaic Dataset which works very well but is there a similar format for vectors?

Mapinfo has hotlinks so we use the index tiles to open files in the frame without having to find the tiles but the best system is the Global Mapper Catalog (gmc) format which works great for vectors and rasters but a search online shows nothing for ESRI (the last mention of this though was in 2005).

Any ideas?


texts:
I would like to add a vector layer to QGIS using the python console (and eventually in a separate python IDE).  I have run the following code from the pyqgis cookbook:

layer = QgsVectorLayer(r'//Users//rwegener//repos//correlations//co_energy_weather_shpfiles//co_energy_weather.shp',
                   'degree_days', "memory")
if not layer.isValid():
    print('layer failed to load')

The layer loads without a problem, however the layer is not valid.

The only other pages I have found on this suggest the user to mess around with slashes so I have also tried the following ways to input the layer:

layer = QgsVectorLayer(r'/Users//rwegener//repos//correlations//co_energy_weather_shpfiles//co_energy_weather.shp',
                   'degree_days', "memory")
layer = QgsVectorLayer(r'/Users/rwegener/repos/correlations/co_energy_weather_shpfiles/co_energy_weather.shp',
                   'degree_days', "memory")
layer = QgsVectorLayer(r'\\Users\\rwegener\\repos\\correlations\\co_energy_weather_shpfiles\\co_energy_weather.shp',
                   'degree_days', "memory")

When I load the file into QGIS with the "add vector layer" button there are no problems.

I have tried also with the cookbook method:

layer = iface.addVectorLayer(r'//Users//rwegener//repos//correlations//co_energy_weather_shpfiles//co_energy_weather.shp',
                   'degree_days', "memory")

only to get the error:

AttributeError: 'NoneType' object has no attribute 'addVectorLayer'


texts: River boundary (bank) extraction from LiDAR

texts: Is there any way to extract a river boundary (bank) from LiDAR data and its products (DEM, intensity image and slope) automatically? Digitizing is very time consuming.


texts:
Using an existing feature, is there a way/method to get the layer which contains that feature?

Example:

# new empty feature which doesn't belong to any layer
new_feat = QgsFeature()

#-------------------------

layer = iface.activeLayer()
features = layer.getFeatures()
# An existing feature which belongs to a layer
feat = next(features) # feat is connected to active layer

# At this point, how to get layer using feature?
# For example:
print(feat.layer()) # or something


texts: How to fill out a field with a matching between an existing field value and a file name in QGIS?

texts:
I'm facing a big challenge. I need to fill out new fields using "Field Calculator" in QGIS. What do I need to load in those fields?: URL's containing the path to a image file (we can use the statement 'base_file_name'). The thing gets complicated when I'm talking about thousands of records (so I cannot copy and paste every single path into the field, I need to make it automatically during a Graphic Modeler processing).

The drop that spilled the glass is that I need that the fields should be filled out with a matching using an existing field in the Shapefile containing a string with a long text (example: "bf5ae62b-a7c9-4682-9c20-4c457de38a9b") and using this code to match it with a file within a folder (each image file name contains this code, for example: bf5ae62b-a7c9-4682-9c20-4c457de38a9b.jpg). I need to find an expression to make it possible and load all the URL at once using this criteria. Thanks in advance if you can help me.


texts: Set encoding when using ogr2ogr to export from Oracle db to MapInfo table?

texts:
I am trying to automate the export of spatial data from a oracle database to mapinfo .tab files. Using this thread I made my script. However I can't see how to set the encoding right.
My data have latin1 encoding. But this is not the case in the output file. I get warnings regarding invalid characters when running the script and Local characters like æøå is not intact in the output.
This is my code:
REM Unsuccessful attempt to change encoding:
SET PGCLIENTENCODING = LATIN1
SET OGR_FORCE_ASCII = NO

C:\OSGeo4W\bin\ogr2ogr.exe -f "MapInfo File" C:\OUTPUT\PATH\test.TAB OCI:"<USERNAME>/<LOGIN>@(DESCRIPTION = (ADDRESS_LIST = (ADDRESS = (PROTOCOL = TCP)(HOST = <HOST>)(PORT = <PORT>)))(CONNECT_DATA = (SID =<SID-NUMBER>))):<TABLES>


texts: Running Linear Regression with multiple Rasters converted to a numpy array in Python


texts: How to calculate interpolated timestamps of linearly referenced points in ArcGIS?

texts:
I have a point feature class with timestamp attribute belongs to bus-gps smartcard data-set. I've used linear referencing tool to snap these points to a bus-route. So there are milepoints in the attribute table of points.
I have another point feature class belongs to bus-stops and i've used linear referencing tool to snap these points to the bus-route. So there are milepoints in the attribute table of bus-stops.

I've merged bus-stop points (having no timestamps) between these gps points (having no timestamps) in the same point feature class.

I need to find the timestamps of bus-stop points between gps points by interpolating according to the distance and time differences (speed) of gps points?

Is there a ready to use tool in ArcGIS toolbox? or can i calculate it in another way?

For example in the following figure, i need to interpolate timestamps of 152, 153, 154 according to the length/speed of 150 and 155.

  Simply what i want to calculate is if timestamp is T and measure is S,

    T152=S152*(T155-T150)/(S155-S150)

      Since there are lots of point features, i need to make it with a tool or script in model builder.


texts:
Returning to QGIS after a long absence.

Probably missing something basic, but I can't click on the Toggle Editing Mode (pen symbol) when I open an Attribute Table in QGIS 2.8.

Any pointers?


texts: How to split single address field into 4 fields (House No., Direction, Name, Suffix)

texts:
I have to do some massaging of our parcel data to make it usable by a program in sheriff helicopters. The program requires one of the following address formats within the fields:

Our addresses are currently in one field:
ex: 1234 W Main St.

Is there a way to automate the splitting of the fields into either of these desired formats?

I can imagine the two field format would be easier by just calling for a split after the numbers, but could also cause a problem for streets such as 1st Ave, etc.


texts: Is it possible to subdivide scale bar segments in QGIS?

texts:
How can I subdivide the scale bar segments in smaller subsegments?

I want to style my scale bar like this:

To make the question clear I edited the posting.


texts:
I recently ran the 10.1 sp1 and my scripts stopped working.  I can't even right click on the scripts and check properties without the following popping up:

InitializeParameters Syntax Error:  Runtime error
Traceback (most recent call last):
File “#ANNOBO2.InitializeParameters.py”, line 1, in module
File “d:\arcgis\desktop10.1\arcpy\arcpy\__init__.py”, line24, in module
From arcpy.toolbox import*
File “d:\arcgis\desktop10.1\arcpy\arcpy\toolbox.py”, line 485, in <module>
From management import ExportTileCache as ExporttileCache_management
ImportError:  cannot import name ExportTileCache

I did a little digging and found all the tools in python script in the management.py.  What was missing was the scripts for the Tile Cache toolset.  Not sure what happened when I installed sp1 but now I just want it fixed.

Any help would be greatly appreciated.

Currently using Python 2.7

I'll look and see if there are any remnants of older python on my system.  If not I'll try to reinstall everything.


texts: Coordinates are not displayed only attribute table

texts: I am new to QGIS and I have uploaded a second set of coordinates as a delimited text layer. I followed all the steps and the coordinates are in the same format as a previously loaded layer, but the new ones do not show up on the map and are shown as layer not supported.


texts:
I am currently using the Map Tip function to display an image when I hover over a point within the canvas. I would like to be able to click on the image to open its folder location. How can I modify the code below to include href and make the image clickable? I would prefer to avoid using the "Actions" function.

%path% = image path

%folder% = folder location

This is the working code I'm using to display the image.

I was trying to use href within the code and I was not successful. <a href="file:///[%folder%]"


texts:
I have a huge point dataset with many different type of fields (text, int and double). I would like to dissolve them by a common ID and for the fields I would like to populate with the most common value (MODE). So far I couldn't find a tool that would have this option.

I am handy with python scripting.


texts: Represent some value along time with QGIS

texts:
I would like to create a dynamic representation of a dataset along the time in QGIS.

I have my data in a layer and a table with a relation from one to many:

TABLE_1      TABLE_2
id_city ---> id_city
city         year
             population

My idea is to make some kind of animation that represent the values of population with a range of colors along the time.

To do that I've created a new layer with the geometry of every city duplicated as many times as years and the data of TABLE 2, but I have some doubts...

Anyone knows if is there any way to do this representation with QGIS?
Is it possible to do it just seting the relation and without a new layer with duplicated geometries?

I have tried the Time manager plugin of QGIS, but as long as i know it shows new geometries along the time, and does not show the same polygon with different values of population (in this case...).


texts:
I have downloaded the QuickMapServices Plugin for QGIS, however when I try to open an OSM layer (or any layer) nothing happens. It appears in the legend, but not in the project workspace (other layers are working fine).
I also cannot access the Contributed Pack, as it says there is an error in the python code.

If anyone else has had this problem, is there a solution?


texts: Road displays up to node, not up to starting point and ending point

texts:
SELECT id FROM table ORDER BY st_distance(the_geom, st_setsrid(st_makepoint($slng,$slat), 4326)) limit 1
; get source and similar for target
; and then apply in dijkstra algo
SELECT pgr_a('SELECT id AS id, source, target,length AS cost, x1, x2, y1, y2 FROM ways',$source,$target,false,false))
SELECT id_0, st_AsText(geom) FROM table w, a a WHERE w.id = (a.pgr_a).id2

and getting linestring but it is not perfect up to staring point and ending point...
it is display up to node not up to point b


texts:
I'm generating elevation contours from the latest ASTER DEMs, which involves processing hundreds of tiles. Unfortunately, some of them have substantial errors, as seen in this colour relief from S01W079:

These patches are not voids. They have "valid" elevations that are just wildly inaccurate. The resulting contours (generated with gdal_contour) look like this:

I need to find these spots, but visually inspecting every single tile for errors is not an option. Is there any method I can use to automatically identify problem tiles? CLI and Pythons solutions preferred, but anything would help.

Right now the only option I see is to write a Python script that picks up strange jumps or drops in elevation across the DEM.


texts: How to open the "File Browse" dialog from PyQGIS plugin code?

texts:
I am in the process of developing a plugin for QGIS using Python. Note that this is my first plugin, so the query could be quite basic...

I have used this Digital Geography link to create the basic plugin, and (quite surprisingly found it easy) created the plugin files.

I am stuck at a very basic concept. Previously, while using Python IDLE, using

QtCore.QObject.connect(self.Input, QtCore.SIGNAL("clicked()"), self.OpenBrowse)

to signal a click and writing a function OpenBrowse

    def OpenBrowse(self):
        self.Input_TB.setText(filename1)

I was able to open the Windows File Browse window on click and set the File location to an adjoining text box.

How do I do this while developing the plugin in the Python file that is created by Plugin Builder.? I am very confused about where should the code be written? Because if I place the 'clicked()' actions, etc. in the Python file generated from the UI file, it is not recognizing the variables as they are not global (obviously).

I am using QGIS 2.4, on a Windows 32-bit system.


texts: Clustering neighboring points inside ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10x.
I have a shapefile of 300 points which I need to make groups, each of which will contain exactly 9 points which are closest neighbors.

I tried the 'near' tool, which have added NEAR_FID and NEAR_DIST to the attribute table but couldn't use it to group 9 neighbors. I also tried 'grouping analysis' tool but this only makes 2 groups overall with variable number of members each time.

Anyone have any better idea?


texts:
I am new in OSM, I have downloaded the data of my map from OSM using the JOSM and i have saved it in a *.xml format.

I'm trying to convert a nodes and a edges (ways) from a OSM map to a graph ( connected graph for example) ...  ... can we do somthing like taht ?? Is there a software to do it ???

someone will say me, try to use the 'GraphView' plugin, but, whaht i need is a specifique software or another way to draw a graph frome xml file for a OSM map ...

I googled a lot, but either I used the wrong keywords or there are no simple solutions on the internet.


texts: Leaflet: IE11 freezes when showing polylines

texts:
I am a newbie in leaflet programming, and I have an error on showing polygons on IE 11.
My page is very simple and works on Chrome and FireFox. I just draw a point and a polyline.
Expected result is shown in the end of the post.All is fine when I just draw a marker but when I add a polyline, on IE, navigation is fine since the screen map doesn't reach the polyline (high zoom level).
As soon as I zoom out to see the polyline, IE11 freezes.
An error is raised in the https://unpkg.com/leaflet@1.6.0/dist/leaflet.js script (5,61705):
[Impossible to get the property "x" of a null or undefined reference]
on function:
// return closest point on segment or distance to that point
function _sqClosestPointOnSegment(p, p1, p2, sqDist) {
var x = p1.x, // Error raises here
    y = p1.y,
    dx = p2.x - x,
    dy = p2.y - y,
    dot = dx * dx + dy * dy,
    t;

if (dot > 0) {
    t = ((p.x - x) * dx + (p.y - y) * dy) / dot;

    if (t > 1) {
        x = p2.x;
        y = p2.y;
    } else if (t > 0) {
        x += dx * t;
        y += dy * t;
    }
}

dx = p.x - x;
dy = p.y - y;

return sqDist ? dx * dx + dy * dy : new Point(x, y);

}
My code is:
function init() {
  var Produit4 = {
    lat: 52,
    lon: 0,
  };
  var ZoomLevel = 7;
  var mymap = L.map("mapid").setView([Produit4.lat, Produit4.lon], ZoomLevel);
  var MainLayer = L.tileLayer(
    "https://api.mapbox.com/styles/v1/{id}/tiles/{z}/{x}/{y}?access_token=sk.eyJ1IjoiYnJ1bm9nYXJuaWVyIiwiYSI6ImNrYmFzam9hdzA5OW0ycnB2Z2JxeXEwYmcifQ.m0w2z_IIm_M2JdO0W7mJBA",
    {
      attribution:
        "Map data &copy; OpenStreetMap contributors, " +
        "CC-BY-SA, " +
        "Imagery © Mapbox",
      maxZoom: 18,
      id: "mapbox/streets-v11",
      tileSize: 512,
      zoomOffset: -1,
    }
  );
  L.Browser.touch = false;
  MainLayer.addTo(mymap);
  var rectangle = L.polyline([
    [46, -16],
    [62.75, -16],
    [62.75, 13],
    [46, 13],
  ]).addTo(mymap);
  var marker = L.marker([51.5, -0.09]).addTo(mymap);
}

Expected result:
"


texts: Making hexgrid in QGIS with cell identification numbers?

texts:
I'm getting started with QGIS as part of my interest in military science and wargaming.  I'm running into a problem adding identifiers (column and row) to the hexcells that I've created with mmqgis.

This image gives an example of what I want to include.


texts: How to link array values stored in a layer to values stored in another layer

texts:
I am stuck for several days trying to get array values link together. I spent hours without getting serious results.
I have a bus stop points table with arrays containing the bus routes numbers serving the stop:
ID transfers
1  (1,2,3)
2  (2,4,A)
3  (1,3,C)
Etc

I have another table with individual values containing the route number and the color
ID    Route     Color
1       1       #CD5C5C
2       2       #FFA500
3       3       #00FF00
4       4       #6A5ACD
5       A       #1E90FF
6       B       #A52A2A
7       C       #008B8B
Etc

I am trying to link them together to get the routes and colors values for each stops as follow:
For example with the stop 3:
Route 1 #CD5C5C to "destination1"
Route 3 #00FF00 to "destination3"
Route C #00FF00 to "destination7"

I have been able to link the colors, but I can't manage getting the route number as well next to the color.
This way I can already get the colors, but I can't link routes for each ones, the best I had was all the routes repeated every time for each color.
array_to_string(

-- Define the variable for transfers Arrays
with_variable('transfers',string_to_array("transfers",','),

array_foreach(
-- Define variables to be used
-- Define the variable for Array index value result
with_variable('index',

-- Define the variable for routes Arrays
with_variable('routes',
aggregate(
layer:= 'Bus_routes',
aggregate:= 'array_agg',
expression:=  "route" ),
-- Find position in list for each route
array_foreach( @transfers, array_find( @routes, @element ) ) ),

-- Get routes color value from array index
with_variable('fill',
aggregate(
layer:= 'Bus_routes',
aggregate:= 'array_agg',
expression:=  "color" ),

array_foreach( @index, array_get( @fill, @element ) ) ) ), @element ) ), ' ' )

EDIT
I partially succeeded using maps arrays. But there's something weird.
If I enter this expression, I won't get the map array, just an empty value.
with_variable('transfers',string_to_array("transfers",','),

with_variable('routes',
aggregate(
layer:= 'Bus_routes',
aggregate:= 'array_agg',
expression:= '''' || "route" || '''' || ':' || '''' || "hex" || '''' ),

map( replace( replace( replace( trim(array_to_string( @routes ) ), ',', ''), '''''', ''':'''), ':', ',') )) )

But if I remove the map command and copy the result from the attributes table and pate it in the code and create a map over it, it works!
For example, this won't work because it gives an empty map result.
with_variable('transfers',string_to_array("transfers",','),

with_variable('routes',
aggregate(
layer:= 'Bus_routes',
aggregate:= 'array_agg',
expression:= '''' || "route" || '''' || ':' || '''' || "hex" || '''' ),

array_foreach( @transfers, map_get(
map ( replace( replace( replace( trim(array_to_string( @routes ) ), ',', ''), '''''', ''':'''), ':', ',') )
@element )  || ' ' || @element) ))

But when I take the values gotten from the replace command in the attributes table and apply a map array over, it will work! That's crazy.
with_variable('transfers',string_to_array("transfers",','),

with_variable('routes',
aggregate(
layer:= 'Bus_routes',
aggregate:= 'array_agg',
expression:= '''' || "route" || '''' || ':' || '''' || "hex" || '''' ),

array_foreach( @transfers, map_get(
map ( '1','#00A440','2','#EA001B','3','#006DB1','4','#FDCA00','5','#00B8E4' ),
@element )  || ' ' || @element) ))

Result: [ '#00A440 1', '#006DB1 3', '#FDCA00 4' ]


texts: How to Export a Shapefile as KML Using R

texts:
All,

I have exported an ESRI shapefile (polygon) to KML using R, but the resulting layer shows only the outer extent of the polygon.
How can I fill in the area with a color?
Any help is greatly appreciated.

The code is as follows:

library(rgdal)
setwd("G:\\GIS_SJR\\GIS Shapefiles\\Boundaries")
polygon <- readOGR(".", "20km_buffer_albers")
polygonWGS <- spTransform(polygon, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
writeOGR(polygonWGS, dsn="polygonWGS.kml", layer="polygonWGS", driver="KML")


texts: I have a point shp file of 445 points around the globe, that i need to create buffer-zone of 10 km around each point, the points shp file is in wgs1984 geographic coordinate and wgs1984 platte carre as projected coordinate system. Every time i create the buffer zone it gives me a HUGE circles above each other that cover the area of whole points together and more!. Does anyone knows why is this happening?


texts:
I imported an OSM file of San Francisco and I get this dataset:

Can any one help me to understand the meaning of these columns:

km/h: is this the limit speed or the average speed of each road?
flags
clazz


texts: Define a action with Python directly (QGIS)

texts:
I'm using QGIS 3.4.3 and I have created some actions for vector layers. I did it interactively and this part is working.

Now I wanted to know if there is a mechanism how I could create my actions directly in Python. I don't want to create these actions interactively in QGIS no more.

I figured out I can handle it with QgsAction. It works for me.

    def my_test(self):
    aLayer = self.iface.activeLayer()

    file1 = self.plugin_dir + "/action1.txt"
    fw = open(file1, 'r')
    txt = fw.read()
    fw.close()

    action = QgsAction (QgsAction.GenericPython, 'test', txt)
    my_scopes = {'Field'}

    actionManager = aLayer.actions()
    i=0
    for a in actionManager.actions():
     if a.name() == action.name():
      i=1
      break;
    if i == 0:
     action.setActionScopes(my_scopes)
     actionManager.addAction(action)

Now I need to know how I can set flag: "Enable only when editable". The option is disabled per default. I tried to use QgsAction with parameter enabledOnlyWhenEditable=True but it doesn't work.  TypeError: 'enabledOnlyWhenEditable' is an unknown keyword argument


texts: Downgrading WRF Hydro GIS Pre-Processing Tools 2.2 to work with ArcGIS Desktop 10.0?

texts:
I have just come across the WRF Hydro GIS Pre-Processing Tools for the first time.  The version I am looking at is 2.2 and its documentation says:

  The processing workflow for creating WRF-Hydro routing grids is
  available to users as an ArcGIS Python Toolbox. Python toolboxes were
  implemented in ArcGIS Desktop version 10.1 as a way to create custom
  geoprocessing tools directly from Python scripts.

Unfortunately, the people wanting to use these tools only have access to ArcGIS Desktop 10.0 which means that they will not "just work" because Python toolboxes did not exist until 10.1.

I am considering trying to downgrade the tools to 10.0 by creating a standard toolbox of Python script tools that calls the same functions that the tools in the Python toolbox do.

Before doing that I am wondering if anyone knows the history of these tools, to be able to tell me whether there was an earlier version of them that were designed to use with ArcGIS Desktop 10.0?  If so, I may be able to try and source that version instead but my Google searches so far have been unsuccessful in finding such a version.

If I need to proceed with the creation of an equivalent standard toolbox for 10.0 then I plan to check that no critical tools/parameters that they use were added at 10.1 and keep my fingers crossed that I have not missed any gotchas.


texts: Is there a way to include a google charts panel alongside cartodb map?

texts: I have successfully included charts in infowindows (e.g. https://github.com/CartoDB/cartodb.js/blob/develop/examples/infowindow_with_graph.html) but I would like to include charts in panels alongside a cartodb map. Is there a way to adjust the code in the link above to achieve this?


texts: I am new to 10.0 and have never used the thiessen polygon tool -- i have a set of points and would like to create a thiessen polygon; what do i need to input (i.e. what type of input (.shp etc)?


texts: Combining gridded polygons at state level within same shapefile using ArcGIS Desktop?

texts:
I have a polygon .SHP containing a field of separate, gridded polygons, which comprise the US and its outliers. I need those separate polygons grouped by state in the same .SHP. I have tried several avenue to join this layer to a state boundary file to no avail.

Ideally I am looking for a manual process.


texts: How to update ArcGIS 10.0 extension for ArcGIS 10.1?

texts:
I have an ArcGIS 10.0 extension which will no longer install for users who have upgraded to ArcGIS 10.1. Are there a defined set of changes that need to be made in order to upgrade the extension, or is it dependent on extension options, such as the libraries that are called?

The error recieved upon install is: "Error 1001. Unable to get installer types in the C:\...extensionname.dll assembly. --> Unable to load one or more of the requested types. Retrieve the LoadedExceptions property for more information."

Any other comments/advice on moving extensions to 10.1?


texts: How to add points/lines/polygons to empty layer in Carto Builder?

texts:
I created an empty map in Carto Builder.  It contained an empty layer.  When I selected that layer in the sidebar, a toolbar for adding points/lines/polygons appeared in the lower left of the map layer.  I used that to add a point.  So far so good.

Then I wanted to add a line.  So I added an empty layer to my map, and selected it in the sidebar.  But the toolbar for adding points/lines/polygons did not appear in the lower left of the map layer.

How do I add a point/line/polygon to a new empty layer in Builder?


texts: Shapefiles intersection with GDAL Python

texts:
I want to compute the area of the intersection for each geometry in two different shapefiles in Python with GDAL/OGR. Here is my sample code:

def ogrIntersection(srcShpPath1, srcShpPath2):
    # import ogr
    driver  = ogr.GetDriverByName("ESRI Shapefile")
    srcShp1 = driver.Open(srcShpPath1)
    srcShp2 = driver.Open(srcShpPath2)
    srcLay1 = srcShp1.GetLayer()
    srcLay2 = srcShp2.GetLayer()
    outArea = []
    for feat1 in srcLay1:
        geom1 = feat1.GetGeometryRef()
        for feat2 in srcLay2:
            geom2 = feat2.GetGeometryRef()
            if geom2.Intersects(geom1):
                outArea.append(geom2.Intersection(geom1).GetArea())
    return(outArea);

ogrIntersection(srcShpPath1, srcShpPath2)

The function returns a single number, while the code iterate over the geometries of the two input layers. What did I do wrong?


texts: Converting KML files for use with Python library Shapely

texts:
I am trying to get a KML-file (wijken.kml) into Shapely. The KML-file is validated against the proper XML schema, so I guess the input is correct.

Routes I have tried:

1) Converting to WKT or WKB format, and reading in with build in functions

Conversion:

ogr2ogr -f CSV wijken.csv wijken.kml -lco GEOMETRY=AS_WKT
ogr2ogr -f SQLite wijken.wkb wijken.kml

In shapely:

from shapely import wkt, wkb
f = open('../kml/wijken.wkb')
wkb.load(f)

Which gives (same for wkt()):

ReadingError: Could not create geometry because of errors while reading input.

As it provides no further info, and the Python code wraps other (C?) libraries, I do not know what is wrong with the format. The CSV-file contains multiple columns, perhaps something is wrong there, but I have not found a comparable online WKT-example to test.

2) Converting to GeoJSON and using the build in Shapely asShape function

ogr2ogr2 -f GeoJSON wijken.json wijken.kml

In Shapely:

import json
from shapely.geometry import asShape
f = open('wijken.json', 'r')
js = json.load(f)
f.close()
asShape(js)

Which gives:

ValueError: Unknown geometry type: featurecollection

This error is the same for a minimal, valid GeoJSON example. Looking at the Shapley code, the problem is that basic GeoJSON types such as "Feature" and "FeatureCollection" are not recognised. It is not clear which route to take to get from GeoJSON features to features Shapely understands.

3) Read in KML with fastkml, which returns Shapely objects

This works, but it seems to lose the properties/ExtendedData values in the KML (or at least does not pass it to the Shapely objects).

I feel that I'm missing something, it can't be so hard to get data into Shapely. Can someone point me towards the most robust/a working solution?


texts:
I am using QGIS to create a frequency map (how frequently does a bus/train run on a route~ every 15mins, 20, 30, 45 or 1 hr) of Greater Cleveland's RTA, from a GTFS fileset (latest version from Transitfeeds). I opened the shape, stops and routes files on QGIS, but cannot extrapolate the individual routes from either files or determine the frequency of the routes.

How can I should proceed?


texts: How to get Coordinate axis (x,y) parallel to orientation in QGIS?

texts:
I am using QGIS 2.14.11. I have a map with drains as you can see on the map. I would like the coordinate (x,y) of each point but my origin point (at the corner just below in red on the map) will be (0,0). The x axis will be parallel at the length and y axis parallel to the width as you can see on the image.

Do you know how that could be done?


texts: Getting features using objectID in ArcGIS Runtime SDK for iOS?

texts:
In ArcGIS Server 10.2, there is a method named lookupFeatureWithObjectId in the class AGSFeatureLayer.h.

However, there is no such method in ArcGIS Runtime SDK for iOS version 100.x.

How can i manipulate such functionality using new version of ArcGIS Runtime SDK for iOS?


texts: How to extend OSM2PO with ferry durations?

texts:
I would like to have the ferry duration [1] taken into account, when transforming OSM files to SQL.

But I'm not sure if I need to write my own waytagresolver, or it is possible to simply add it to the osm_meta tag by adding a parameter or by config? I'm hoping for the latter, as this of cause is the simples by far :)

[1] http://wiki.openstreetmap.org/wiki/Tag:route%3Dferry


texts: Make linear regression in multiple raster layers

texts:
I have a question about multiple raster layer regression using R or Arcgis.
I would like to predict biomass using some Vegetation Index - VI (NDVI, EVI, SAVI, OSAVI...). and generate R² and RMSE maps for my study area.

At first moment, I cut the VI images and stacked them with biomass, because I have for some areas biomass field plot. So, I create a stack with biomass field and VI's. Then I extracted the value (Biomass and VI's for each pixel where I have biomass value) and generated a model to predict biomass (I have used R and regsubset function to do it). Doing this a got, for example (Biomass = B0 + B1*NDVI + B2*SAVI) I could apply this model in my VI's images to extrapolate, but I will have only biomass maps.

I think to get R², RMSE maps, I must do a different process, someone could help me?


texts:
I'm looking at replacing our corporate GIS, this is the first time that I have had the responsibility of such a task despite my years of experience with the software as an end user.

We will be looking at an ESRI platform for our desktop administrative functions but in terms of storage, database structure etc I'm beginning to realise what a minefield this is and hope that I can draw on the knowledge and experience of you good people to perhaps provide some guidance in what to look for perhaps highlight any potential pitfalls I may stumble across along the way.

My current questions are :-

FGDB vs SDE (which will mean having to puchase ArcGIS for Server).  How can I best justify the investment of server just for the purposes of data storage (we will be using a different vendors product for serving up our GIS data).
POSTGRESQL/POSTGIS vs ArcGIS for Server - is postgresql/postgis a direct alternative to ArcGIS for Server/SDE etc or are the two completely different...and if so how ?
Is it easy and straightforward to connect to POSTGRESQL/POSTGIS from within ArcGIS Desktop and edit data or does this require some other 3rd party software such as FME or GISSQUIRREL etc ?

My preference is to have ArcGIS for Desktop and ArcGIS for Server from a continuity perspective.  There is also the other issue that we are looking at integrating with our other corporate systems, so ideally we only want to do this once and my concern is that if we have different bits of software for different purposes that we are going to end up with more integrations that we perhaps need.

I just want to try and avoid extra steps in processes here, from data creation, to maintenance to publishing etc etc and whilst Dekho would have been the logical step given the fact we will be using ArcGIS Desktop and potentially Server, we felt the end user (GIS desktop viewer and field viewer) requirements we best fulfilled elsewhere so we are happy to take the hit on the extra steps involved in serving our data up this way...which in all honesty isn't that big a deal as most of it is automated anyway.


texts:
I have an extern geojson file (points) and I would like to display a circle for each point coordinates and get a bigger circle on a mouse click, with a popup.

I have no problem to load the geojson and to view the circle markers, but I can't get both bigger cricle and popup.

Here is part of my code :

function clickFeature(e) {
                    if(clickedFeature) {clickedFeature.setRadius(3)}
                    var layer = e.target;
                    e.target.setRadius(5);
                    clickedFeature = e.target;
                    info.update(layer.feature.properties)
                    }
function onEachFeature(feature, layer) {
                                    layer.bindPopup(feature.properties.PROPERTY).addTo(MyLayer);
                                    layer.on({click: clickFeature,})
                                    }

MyGeojson = new L.geoJson(mygeojson, {
    onEachFeature: onEachFeature,
    pointToLayer: function (feature, latlng) {
    return L.circleMarker(latlng, {
                radius: 3,
                fillColor: getColor(feature.properties.TYPE),
                color: "#242121",
                opacity: 1,
                weight: 1.5,
                fillOpacity: 1
            });
        }
    }).addTo(map);

First I had this same kind of functions and code, but using PNG icon instead of circleMarkers, and it worked well, both bigger icon and popup on mouse click.

Any idea what is wrong ?


texts: Which provider should I use with `qgis.core.QgsVectorLayer` method at `PostgreSQL database with PostGIS extension` storage?

texts:
I use PyQGIS 3.6.1 Noosa to load layer from my database storage PostgreSQL database with PostGIS extension using qgis.core.QgsVectorLayer method. But layer is invalid somehow. Performing same operation with QGIS interface - successfully creates the layer (see picture).

How to create a layer with qgis.core.QgsVectorLayer method?

Here is part of my code:

def load_layer_from_DB(aSchema, aTable, aGeometryColumn, host, port, dbname, user, password):
    uri = qgis.core.QgsDataSourceUri()
    uri.setConnection(host, port, dbname, user, password)
    uri.setDataSource(aSchema, aTable, aGeometryColumn)
    layer = qgis.core.QgsVectorLayer(uri.uri(), 'some_name', "postgis") # i have tried `postgres`, `postgresql` but no luck

# connect stuff
dbname = 'db_visdata'
host = 'ххх.ххх.ххх.ххх'
port = 'хххх'
user = 'хххххххххх'
password = 'хххххххххх'

# path in database
aSchema = "test_data"
aTable = 'test_satellite_boarders'
aGeometryColumn = 'MultiPolygon'

layer = load_layer_from_DB(aSchema, aTable, aGeometryColumn, host, port, dbname, user, password)

print('layer.isValid=', layer.isValid())
print('layer.type=', layer.type())
print('layer.name=', layer.name())
print('layer.featureCount=', layer.featureCount())
return layer

Output:

layer.isValid= False
layer.type= 0
layer.name= some_name
layer.featureCount= -1


texts: Moving a vector layer 20m to the west on QGIS

texts: I have a vector layer of power lines for an area, the shapefile that I was given is actually around 20m to the west of where I know the power lines are actually supposed to be, as they are supposed to coincide with the main road and they do not. Is there an easy way to move this vector layer slightly to the east? I am using QGIS.


texts:
I have a problem. MapServer displays raster on top other layers. Is there any way to put raster beneath the map?
MAP
  IMAGETYPE      PNG
  EXTENT         466914.4981 5983965.6134 469588.2900 5986580.8550
   SIZE           300 300
  SHAPEPATH      "../data"
  IMAGECOLOR     255 255 255

     # Start of LAYER DEFINITIONS ---------------------------------------------
  LAYER # States polygon layer begins here
    NAME         "Pelke_rastras"
    DATA         "raster\Pelke_rastras4.tif"
    STATUS       ON
    TYPE         RASTER

   CLASS
      NAME       "Krakinis 1995"

      END
    END
  END # States polygon layer ends here

  # End of LAYER DEFINITIONS -------------------------------

  # Start of LAYER DEFINITIONS ---------------------------------------------
  LAYER # States polygon layer begins here
    NAME         Pelkinis_miskas
    DATA         Pelkinis_miskas1
    STATUS       OFF
    TYPE         POLYGON

   CLASS
      NAME       "Krakinis 1995"

      STYLE
        COLOR        102 204 0
        OUTLINECOLOR 76 153 0
      END
    END
  END

 Ežeras "Krakinis" 1995 m.

            <img border="1" src="http://localhost/cgi-bin/mapserv.exe?map=/ms4w/apps/tutorial/htdocs/3LAB.map&layer=Pelke_rastras&layer=Pelkinis_miskas&layer=Pelke&layer=Ezeras&mode=map"


texts: Esri Roads and Highways Extension on Linux ArcGIS Server

texts:
Is it possible to add Esri Roads and Highways Extension on ArcGIS Server running on Linux?
The products list only includes windows version.


texts:
I have a requirement to convert .KML file to GeoJSON or JSON. Can you please recommend open source command line tools on GNU/Linux(2.6.32-358.14.1.el6.x86_64)? I have installed gdal, when I run..

ogr2ogr -f "GeoJSON" geojson.json Areas.kml

I get this error :-

  Unable to open datasource `Areas.kml' with the following drivers.   ->
  ESRI Shapefile   -> MapInfo File   -> UK .NTF   -> SDTS   -> TIGER  -> S57   -> DGN   -> VRT   -> REC   -> Memory   -> BNA   -> CSV   -> GML   -> GPX   -> KML   -> GeoJSON   -> GMT   -> AVCBin


texts: Layer not loading correctly after table join

texts:
I am attempting to perform a table join from a plug-in/outside of the built-in Python console. This is QGIS 3.16.1 and Python 3.7. I want to join two vector point layers, which have all points in exactly the same places. This is my code for the join:
# Create layers for output files
vectorpoint_filled_one = QgsVectorLayer(outputfile, '', 'ogr')
vectorpoint_filled_two = QgsVectorLayer(outputfile_two, 'final', 'ogr')

# Join tables
join_info = QgsVectorLayerJoinInfo()
join_info.setJoinLayer(vectorpoint_filled_one)
join_info.setJoinFieldName('field1')
join_info.setTargetFieldName('field1')
join_info.setUsingMemoryCache(True)
join_info.setJoinFieldNamesBlockList(['geom_X', 'geom_Y'])
vectorpoint_filled_two.addJoin(join_info)

# add to map canvas
QgsProject.instance().addMapLayer(vectorpoint_filled_two) #loads version from before table join

For some reason, the loaded map is vectorpoint_filled_two before the table join. However, when I also load vectorpoint_filled_one, the vectorpoint_filled_two is loaded (correctly) after the table join:
# Create layers for output files
vectorpoint_filled_one = QgsVectorLayer(outputfile, '', 'ogr')
vectorpoint_filled_two = QgsVectorLayer(outputfile_two, 'final', 'ogr')

# Join tables - add if statement to skip for no vector layers or no raster layers
join_info = QgsVectorLayerJoinInfo()
join_info.setJoinLayer(vectorpoint_filled_one)
join_info.setJoinFieldName('field1')
join_info.setTargetFieldName('field1')
join_info.setUsingMemoryCache(True)
join_info.setJoinFieldNamesBlockList(['geom_X', 'geom_Y'])
vectorpoint_filled_two.addJoin(join_info)

# add to map canvas
QgsProject.instance().addMapLayer(vectorpoint_filled_one)
QgsProject.instance().addMapLayer(vectorpoint_filled_two) #now loads correctly

This works, but I don't want to load vectorpoint_filled_one, as it is otherwise useless. I have tried various manners of reloading and updating vectorpoint_filled_two after the table join (startEditing() + commitChanges(),  reload()), but it keeps loading the older version. I have also tried removing vectorpoint_filled_one after vectorpoint_filled_two is added to the map, but this will also result in the join missing from the loaded layer. Since the second version does work and addJoin returns True, I know the table join and previous code creating the layers are correct.
How do I load the correct version of vectorpoint_filled_two without also loading the other layer?


texts:
I have a shapefile of polygons. What I want to do is take some areas and discard others. For example, I have a map of USA (shapefile). I want to take California only and make another shapefile.

How do I do it?


texts: Determining UTM zone from GeoTIFF in QGIS

texts:
I'm very new to QGIS, and I have a sample GeoTIFF (and corresponding TFW) file(s) that I'm trying to determine what the longitude and latitude coordinates for each of the 4 corners of the image are - coordinates of 2 opposite corners of the image would also be okay.

I'm able to import the GeoTIFF as a raster layer in QGIS, but I noticed that all the coordinates that are displayed are in UTM format. In the metadata, I see the layer extent, which I was told are the UTM coordinates of the bottom-left and top-right of the corners of the image (square-shaped), which I can work with. However, I have to determine what the UTM zone is so I can use that along with the UTM coordinates in the metadata to convert the points to latitude and longitude.

This seems like a fairly simple task, but I'm unable to find out what the zone is of the image to move forward with the UTM to long/lat conversion.


texts: How to conserve layer blending after exporting to SVG and PDF in QGIS?

texts:
When exporting a map that make use of label blending in QGIS 3, only when exporting as PNG is the blending kept. If I try exporting it to SVG or PDF, the blending is lost.
Is there any way to export to SVG or PDF while conserving the blending mode of layers that make use of it? The image below illustrates the actual (left) versus the expected (right) results.

The only fixes I have found so far are either to manually add the blending afterwards in Inkscape (but it's time-consuming and the text is also blended), or to purposefully overlap two layers, one with the shadow background and one with the text, so that only the shadow background layer can be blended in Inkscape (but I haven't managed to do so in QGIS).


texts: Is there a text box parameter type for Python Toolboxes?

texts:
The Calculate Field tool has a text box for code.

Is there a parameter type for Python Toolboxes for this type? There's String but it's just one row of text.

    in_str = arcpy.Parameter(
        name='test',
        displayName='Test',
        datatype='String',
        direction='Input',
        parameterType='Required')