texts: How to generate smooth contours from cropped raster?

texts:
I’m trying to create contours from a raster layer in QGIS. On my first attempt, I cropped my raster layer in Photoshop and reimported it, but it lost coordinate data I think. I made the contours successfully and they were beautiful, but mirrored upwards from the raster layer.

The scale is way out. I discovered this when I imported some vector data. A distance that should be 0.04 degrees is now 30 degrees. This probably is why I’m having my next problem:

I eventually figured out a method for cropping the raster layer in QGIS so it’s at the right scale and location. However, when I make contours from this raster they appear square and pixellated, like they follow the pixel boundaries of the raster layer.

The pixel data for the two raters I used looks the same, so why is this one at the smaller scale now looking so different? Any idea how I can get the nice smooth contours that I got on my first attempt (albeit at the wrong scale)?

I’m very new to this, so I may be missing obvious things! I’d appreciate any help.


texts: Accuracy of Google Maps latitude and longitude?

texts:
I am making a program that uses GPS location preferably from smart phone, laptop works too and compares it to a building, the points you in its general direction. I am making this web application specifically for San Jose State Campus for fun. I have read the other posts about the correlation between decimal places and accuracy and a .11m error will be sufficient. Google maps gives me coordinates to the 6th decimal place, I was wondering if they do indeed follow the standard I have read about here. Here is an example of one of the buildings I will index.

https://maps.google.com/maps?q=%2B38%C2%B0+34%27+24.00%22,+-109%C2%B0+32%27+57.00%22&um=1&ie=UTF-8&sa=N&tab=wl

I will only have an arrow point in the direction of the building and rotate it as the lat/long of the user changes and maybe display the approximate distance from the building. I have yet to learn how to properly use google api to draw the users path and etc and that is why I am asking.


texts: Is it possible to print the output from ogr2ogr directly to stdout?

texts:
I know I can convert a Shapefile to GeoJSON as the following:
ogr2ogr -f GeoJSON outputFile.json inputFile.shp

Instead of writing the result of this conversion to the file outputFile.json, I'd like to print this result directly on stdout, without having to write a file on the system. Since I'm inside a bash script I know I can find alternative solutions as the following one:
tmpFile=$(mktemp).json
ogr2ogr -f GeoJSON $tmpFile inputFile.shp
cat $tmpFile

However, I'd like to know if ogr2ogr supports printing to stdout instead of writing on files, just so I can do it in one line of code. I've searched in the man ogr2ogr page but it doesn't even contain the word stdout there. Is it possible to print the output of ogr2ogr directly to stdout?


texts: ArcPy command fails (pathname becomes corrupted) based on the starting letter of input file?

texts:
I'm using the Python window in ArcMap 10.3.1 to conduct Maximum Likelihood classifications (MLClassify, Spatial Analyst extension) on multiple images. I've done this with no issues before, but now when I run the code, the pathnames for some (but not all) of the inputs seem to be altered and cause errors.

Here is sample code I'm running in the Python window that defines my input files and runs the classify tool:

import arcpy
from arcpy import env
from arcpy.sa import *

aug1951_clip_BW = "C:\Project\GIS\Grayscale\GrayAllYears.gdb\aug1951_clip_BW"                                #input raster
aug1951_clip_bw_gsg = "C:\Project\GIS\Grayscale\aug1951_clip_bw.gsg"                                         #signature file, used in classifying
Class_aug1951_clip_BW = "C:\Project\GIS\Grayscale\GrayAllYearsNewGSGOutput.gdb\Class_aug1951_clip_BW"        #name of output classified image
ConClass_aug1951_clip_BW = "C:\Project\GIS\Grayscale\GrayAllYearsNewGSGOutput.gdb\ConClass_aug1951_clip_BW"  #name of optional confidence image output

arcpy.gp.MLClassify_sa(aug1951_clip_BW, aug1951_clip_bw_gsg, Class_aug1951_clip_BW, "0.0", "EQUAL", "", ConClass_aug1951_clip_BW)

The MLClassify command then fails to execute--interestingly, in the error message that appears, the "\a" in the signature file (.gsg) pathname seems to be replaced with a black dot/bullet point. My guess was that this altered pathname is what renders the file in the argument unrecognizable and causes the following  error:

  ExecuteError: Failed to execute. Parameters are not valid. ERROR
  000865: Input signature file:
  C:\Project\GIS\Grayscale•ug1951_clip_bw.gsg does not exist. Failed to
  execute (MLClassify).

However, this doesn't happen with my signature files that do not begin with the letter "A"!
Below, an example of a classification that ran with no errors, identical in all ways besides the input filenames:

sep1979_clip1_BW = "C:\Project\GIS\Grayscale\GrayAllYears.gdb\sep1979_clip1_BW"
sep1979_clip1_BW_gsg = "C:\Project\GIS\Grayscale\sep1979_clip1_BW.gsg"
Class_sep1979_clip1_BW = "C:\Project\GIS\Grayscale\GrayAllYearsNewGSGOutput.gdb\Class_sep1979_clip1_BW"
ConClass_sep1979_clip1_BW = "C:\Project\GIS\Grayscale\GrayAllYearsNewGSGOutput.gdb\ConClass_sep1979_clip1_BW"

arcpy.gp.MLClassify_sa(sep1979_clip1_BW, sep1979_clip1_BW_gsg, Class_sep1979_clip1_BW, "0.0", "EQUAL", "", ConClass_sep1979_clip1_BW)

Why is this bullet point inserted in the middle of my signature file pathname/how can I stop this from occurring, and how can I execute MLClassify for these images without having to rename my signature files that begin with "A"? This is the only trend I find in common among the images that have failed to run.


texts: Plotting a value above the line, at a certain distance from begin and end point

texts:
I am using QGIS 3.12

I am trying to find a way to plot a value above the line, at a certain distance from point 1 to the right, and at a certain distance from point 2 to the left ?

The value is coming from a attribute table where I don't want to change or add columns.

So far I used rule-based labelling, placement, data refined....but I'am not there yet... The values are positioned at the point, both begin value and end value, are now presented at the same location.... but needs to presented above next to each other.

How can I do this?

I need to make the label in XML and SLD-format


texts:
I have read the gdal documentation for gdal_rasterize. As well as reading the original forum that established -tap (target aligned pixels). I understand -tap  is used to make a standard grid (maybe like ArcGIS snap raster) but I do not understand the usage.
When I use gdal_rasterize and I enter -tap what is -tap actually doing to the ouput raster?
I do not seem to be able to point -tap to another raster to use this grid and I do not seem to be able to enter origin coordinates. Gdal_rasterize appears to want me to enter only -tap with no following argument. My assumption is it aligns with 0,0 but this is a guess.


texts:
I have an elevation-raster in .hgt format with CRS EPSG:4326. My aim is to calculate aspect and slope which - as far as I understand - are difficult to obtain when using geographic CRS. Projected CRS are much better suited as map units are not degrees but meters (or the like).

So I tried and transformed the layer to EPSG:5677. However, the transformation produces very unpleasant artifacts (see pictures). The angles which were produced by the transformation have the potential to distort further analyses to a significant extent, I fear.

Can anyone suggest a procedure to minimise these unwished effects?

before transformation:

after transformation:

close before transformation:

close after transformation:


texts:
I am trying to write a script (python) to open an external exe file on the server once I click on a feature in ArcGIS. The script must have 2 parameters in namely FacilityId and UserID. The exe file will take the parameters and write it to an SQL table where another function triggers from.

This is the python script I have thus far... (Please forgive, me I have not clue what I am doing):

import webbrowser
def OpenLink ( [FacilityId] ):
  path = "W://QGIS//PREMISQGIS.exe"
  webbrowser.open(path)
  return

The above script at least opens the exe and asks if I want to run it - as soon as I do, it says it doesn't have parameters?


texts:
I want to calculate the values for each feature using python.

First I added a new field and it worked excellently (Adding field and calculating expression with PyQGIS?).

Now I want to calculate an Expression

Here is my code

from PyQt4.QtCore import QVariant
from qgis.core import QgsField, QgsExpression, QgsFeature

#add new field
myField = QgsField ('eur1_km2', QVariant.Double)
layer.startEditing()
layer.dataProvider().addAttributes([myField])
layer.updateFields()
idx = layer.dataProvider().fieldNameIndex('eur1_km2')

#calculating values for myField
e = QgsExpression( '[SHAPE_Area]/1000000*6' )

for f in layer.getFeatures():
    f[idx] = e
    layer.updateFeature( f )

layer.commitChanges()

The "SHAPE_Area" is a Field within the attribute table with different values. It is in square metres.

How can I call the values for "SHAPE_Area", so it can be calculated?

I'm using QGIS 3.0 for this dataProvider() is necessary, otherwise an error occurs.

idx = layer.fieldNameIndex('eur1_km2')
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
AttributeError: 'QgsVectorLayer' object has no attribute 'fieldNameIndex'

The next is, that in QGIS 3.0 pendingFields() does not longer exists. An error occurs

e.prepare(layer.pendingFields())
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
AttributeError: 'QgsVectorLayer' object has no attribute 'pendingFields'

Last point another error arise for the last part

for f in layer.getFeatures():
    f[idx] = e.evaluate( f )
    layer.updateFeature( f )
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 2, in <module>
TypeError: QgsExpression.evaluate(): arguments did not match any overloaded call:
  overload 1: too many arguments
  overload 2: argument 1 has unexpected type 'QgsFeature'

So, is there an other method for calling the attribute values?


texts: Artifacts in NetCDFs served to Leaflet from Geoserver over WMS when zooming

texts:
I'm serving a set of NetCDF layers to a Leaflet map over WMS. The netCDFs are native EPSG:4326, and Leaflet requests them in EPSG:3857. They look great zoomed out:

But if I zoom in, I start seeing artifacts and gaps in the tiles that move as I continue to zoom in:

EDIT: some more details on the artifacts. When they start appearing at zoom 4, they occur in the 'left-most tiles' (whichever ones have a lower-left X value of about -20 000 km) and in the tiles directly to the left of X = 0 (but in the top-left this time). For example, this is over Scotland:

The width has been variable with zoom, but they generally get taller until there's just a continuous strip at higher zooms.

This image suggests that there's some problems going across the entire X < 0 space even though artifacts are only visible on the left-most tile (aligned with its right edge) and the right-most tile (aligned with its left edge):

I note that the gaps between tiles appear transparent, just as my regions with no data (NaN) do. In other regions, I can see the tile boundaries shifting on both sides (so no gap appears) when I zoom in.

The artifacts appear in regions that otherwise wouldn't have any data, and they appear as whatever the lowest-valued colour in my style is, regardless of the style I use.

This behaviour appears in all of the NetCDFs I've tested (a set of about 30 files), which are all produced through the same process. However, if I use the style preview when editing a style (which looks like a different projection to my eye), I don't see any artifacts:

This makes me wonder if some sort of error is occurring in the projection to EPSG:3857, either on Geoserver's part or on Leaflet's in requesting the tiles.

I note that the artifacts also appear in the requested tiles when checking the individual requests; this isn't a problem with their rendering in the browser. I've tested this both on a local installation of Geoserver 2.12.1 on macOS and on a Tomcat 7 installation of Geoserver 2.12.2 on CentOS.

EDIT: forcing Geoserver to serve EPSG:4326 by including srs: 'EPSG4326' in my layer request also seems to fix this (which could work for me, especially as GBIF hosts basemaps in EPSG:4326):


texts: Error reprojecting point cloud with liblas 1.8: "latitude or longitude exceeded limits"?

texts:
I was previously reprojecting a point cloud in .las format using liblas 1.7.0 (with GeoTIFF 1.4.0 GDAL 1.10.1) as follows:

las2las /tmp/input.las --a_srs EPSG:4326 --t_srs EPSG:3857 --scale 0.001 0.001 0.001 --offset -13704940,4639778,0 --output /tmp/new.las

However, after upgrading to liblas 1.8.0 (with GeoTIFF 1.4.1 GDAL 1.11.3) I now get the following error which I believe indicates the scale factor and offset combination are not enough to cover the data range:

  ERROR 1: latitude or longitude exceeded limits
  error: Could not project point for ReprojectionTransform::latitude or longitude exceeded limits

I have played about with the offset and scale factor and can get it to work only if I set the scale factors to 0.01 and remove the offset, e.g.:

las2las /tmp/input.las --a_srs EPSG:4326 --t_srs EPSG:3857 --scale 0.01 0.01 0.01 --output /tmp/next.las

Why am I getting this error with the newer version of liblas, and what should my reprojection command be?

For your information, here is the lasinfo for the input point cloud:

  Version:                     1.2
  Source ID:                   0
  Reserved:                    0
  Project ID/GUID:             '00000000-0000-0000-0000-000000000000'
  System ID:                   'libLAS Python API'
  Generating Software:         'libLAS 1.8.0'
  File Creation Day/Year:      149/2018
  Header Byte Size             227
  Data Offset:                 738
  Header Padding:              0
  Number Var. Length Records:  4
  Point Data Format:           2
  Number of Point Records:     20
  Compressed:                  False
  Number of Points by Return:  20 0 0 0 0
  Scale Factor X Y Z:          0.00000010000000 0.00000010000000 0.01000000000000
  Offset X Y Z:                -123.1140000 38.4250000 -100.00
  Min X Y Z:                   -123.1139429 38.4257008 -2.43
  Max X Y Z:                   -123.1130518 38.4263682 18.19
  Spatial Reference:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]


texts: postgis distance calculation between points with lat, long coordinates

texts:
Team,

I am new to postgis, but a postgres user for quite a while.  I would like to compare

the geospatial capabilities between postgresql and mongodb.

I have the lat.lon, zip code of US locations;

I ran across this blog entry at http://unserializableone.blogspot.com/2007/02/using-postgis-to-find-points-of.html

Also, I see this the function ST_DWithin

SELECT * FROM geotable
WHERE ST_DWithin(geocolumn, 'POINT(1000 1000)', 100.0);

Please advise,

I see that the lat/lon must be converted into points; but I am having issues with the best mechanism to perform indexed distance calculations.

Update..

I have made progress and loaded a spatial aware table .. see below

create table public.disttest
(

    zipcode varchar(5),
    state varchar(2),
    city varchar(55),
    long numeric,
    lat numeric

);

SELECT AddGeometryColumn( 'public', 'disttest', 'geom', 32661, 'POINT', 2 );

COPY public.disttest FROM 'c:\temp\zips.txt'

select * from public.disttest limit 10;

UPDATE public.disttest SET geom = ST_Transform(ST_PointFromText('POINT(' || long || ' ' || lat || ')',4269),32661) ;

Now the question is how to use the  below type query for a 10 mile radiuus

SELECT * FROM public.disttest WHERE ST_DWithin(geom, geom, 16093);


texts: Retrieve available PostGIS connections in PyQGIS

texts:
Can I retrieve the available connections to PostGIS databases in PyQGIS? I would like to provide a list of available db-connections, and subsequently a list of tables within the ui of my plugin.

I checked the cookbook but cannot find a way to get further with this.


texts: HTML Dom elements as marker / feature in Openlayers?

texts: I want to set HTML Dom element as marker in openlayers, instead of GIFs,PNGs. Is it possible to do so?


texts: How to colour polygons with a completely transparent fill using MapServer?

texts:
Is it possible to fill polygons so that they are completely transparent in MapServer? I want to have a black background that is slightly opaque, and wherever there is a polygon I want that to be completely clear with no fill. I set my background in the following way:

   IMAGECOLOR     0 0 0
   TRANSPARENT OFF

I've tried creating a class that is white to fill the polygon, and experimented with opacity using TRANSPARENCY with values from 0 - 100 but it does not result in a clear polygon:

CLASS
      NAME "light"
      EXPRESSION ([max_cost] == 2000)
      STYLE
        OUTLINECOLOR 255 255 255
        COLOR 255 255 255
        SIZE       1
        WIDTH      2.0
      END
END

This is results in the following type of image:

Instead of the polygon having a white fill; I want it to have a fill with no color. I found one suggestion online to define the color like this:

COLOR -1 -1 -1

but it also did not work. In this part of the documentation, the following is described:

  “PALETTE=/path/to/palette.txt” is used to define the absolute path
  where palette colors can be found. This file must contain 256 entries
  of r,g,b triplets for RGB imagemodes, or r,g,b,a quadruplets for RGBA
  imagemodes. The expected format is one triplet (or quadruplet) per
  line, each value separated by commas, and each triplet/quadruplet on a
  single line. If you want to use transparency with a palette, it is
  important to have these two colors in the palette file: 0,0,0,0 and
  255,255,255,255.

After reading this, I also tried using COLOR 0 0 0 0 but this throws an error:

  Unknown identifier. Parsing error near (0):(line 71)

I don't understand how I refer to a define color in PALETTE=/path/to/palette.txt but this seems like one possible way to achieve this.


texts: Using ExportToJPEG to plot multiple raster files with arcpy

texts:
I am trying to use arcpy.mapping.ExportToJPEG to convert some raster files to jpeg format. I am new to this module. I am using the following script:
mxd = arcpy.mapping.MapDocument(r"E:/analysis/maps/new_cs_rivnat/c1_s1/53_cities.mxd")
arcpy.mapping.ExportToJPEG(mxd, r"E:/analysis/maps/new_cs_rivnat/c1_s1/test.jpg")

In the 53_cities.mxd, there are 12 raster files and I want to convert all these, How can I do that? This script only converts the first raster.


texts: Calculating NDVI values for specific zones using zonal statistics using QGIS

texts:
I am currently working with multispectral data from a MicaSense RedEdge M on a plot size study in wheat research plots. I want to conduct zonal statistics to calculate NDVI for individual plots.

I found a YouTube video from MicaSense entitled "Calculating zonal statistics for RedEdge-M data using QGIS" and the workflow helped me understand some things. However it did not go into detail on how to obtain the NDVI values. I have created an NDVI raster layer, by using the raster calculator and of course the formula for NDVI and respected bands. I created a grid layer, placed the polygons over the desired area, and conducted zonal statistics for each polygon.

I do not understand the statistical output in the attribute table neither how to convert the data to an NDVI index. Do you have any information on that?


texts: Uniform Annotation Font Size in Multiple Data Frames

texts:
I'm making a map that has five inset maps that display a smaller map scale where point features are more dense.  I have exported my labels to annotations in the geodatabase so that I could fit the labels that maplex wouldn't draw and move things around for the final touches.  I am using the same annotation layer across all the insets so my GDB isn't cluttered with multiple annotation classes.

The problem is that because these data frames use different map scales the annotations are drawing at different sizes because it is set to a single reference scale.  I can't set it to 0 because it is in a geodatabase.  Using Symbol Substitution changes the font size, but it doesn't scale the font to the same sizes (10pt in one data frame displays at a different size than 10pt in another data frame).

My question is:
How do you display the same geodatabase annotation class in multiple data frames, with diffent scales, at the same page layout font size?


texts: DataB and DataB-copy worked on simultaniously - how merge DataB-copy into DataB

texts:
ArcGis 10 - EDITED (sorry for not being more specific, quite new to the arcgis world)

We need to enter A LOT of Map data (Sewage system) into a Database in a relatively short time. It involves an area of about 14km^2 with pipes, manholes (couple hundred) and relvant buildings. All of which have numerous attributes. These are cataloged in 30 some feature classes. Data entered via scans that are georeferenced and then copied by hand via ArcMap. I assume that they're ESRI-feature classes [ArcCatalog says "Personal-Geodatabase Feature-Class"].
The database (.mdb file) is about 130mb.

Following idea:
Database "A" is main Database. A copy is made "A-copy".
With two liscenses two people can work on A simultaneously (different areas).

After a certain time, we want to import all the entered data from "A-copy" back into "A".
How would one best attempt to do this this without overwriting the changes made in "A" during this time?

(the entered data is as good as idependant from the data already entered; It would be possible to wipe all tables from "A-copy" before starting)

Both users would ideally be accesing and entering data into the same tables. [Maybe it would be possible to limit that to a certain extent, not quite sure how topo-lines owuld react if e.g. one user placed all the manholes in one area and then on the next day a second person would connect them in the copy-db and then transfer the pipe feature-class back into the original.]
From what I can tell sde is not in use; we have no possibility of multi-user simultaneous access to this particular database. Is there a way to qickly check this?
I didn't set up the system and am not versed in doing that either.

This wouldn't have to happen on a daily basis, so a multistep process would be still be feasable.

Thanks a million


texts: PostgreSQL query selection in MapServer is not working?

texts:
When I make a selection in the query in Mapserver I got a blank page.

The query I use is:
DATA "geom from (select * from bagmutlijst_201210.vo_gepland_gemeente where gemcode = '0141') as selection using unique gemcode using srid=28992"

The data is correct because if you use the following query (without a selection) I got a correct map.
DATA 'geom from bagmutlijst_201210.vo_gepland_gemeente using unique gemcode using srid=28992'

Why is the selection not working?


texts: Trigger custom geoprocessing tool outside ModelBuilder?

texts:
I have created a model in a custom toolbox with ModelBuilder in ArcMap 10.1.

Is it possible to trigger the model outside ModelBuilder and ArcMap?

For example, use a Python script to trigger the model to operate?


texts:
Each time I run

gdalwarp -te 212051.9319179999874905 214051.9319179999874905 183962.2946890000021085 Quada.tif Clipped.tif

I get the following error:

FAILURE: No target filename specified.

Any ideas as to why?


texts: Adding Styles like LineColor, Fill Color to SimpleFeaturePolygons while exporting them into a kml file

texts:
I am trying to export features(polygons) in to a KML file. I am able to get KML file with the given features. Here all lines are coming in a white color. How can i give different colors to different polygons. Please suggest me.

Here my code

    SimpleFeatureTypeBuilder builder = new SimpleFeatureTypeBuilder();
    builder.setName(plotName);
    builder.setCRS(DefaultGeographicCRS.WGS84);
    builder.add(sPolygon, Polygon.class);
    builder.length(32).add("AreaLocation", String.class);
    final SimpleFeatureType featureType = builder.buildFeatureType();

getting Feature Collection from featureList

    SimpleFeatureCollection collection = new ListFeatureCollection(builder.getFeatureType(), featureList); //featureList list of all polygons

exporting as a kml

            OutputStream stream = new FileOutputStream("my.kml");
            Encoder encoder = new Encoder(new KMLConfiguration());
            encoder.setIndenting(true);
            encoder.encode(collection, KML.kml, stream);


texts:
I have some MultiPolygon geometries in SRID 4258 and am trying to transform them into SRID 4326. When I run the following command, however, the output still seems to be in SRID 4258.

Command:

SELECT ST_AsText(ST_Transform(ST_Centroid(geom), 4326))
FROM   geoms_to_convert
LIMIT  1

Ouput:

POINT(3819288.74433993 3132617.8860223)

I am converting to a point to see if the transform works, and the end goal is to have the point return in lat and long.

I don't see what I am doing wrong and am relatively new to PostGIS. I have confirmed that both projections are present in the spatial_ref_sys table and am using PostGIS version 2.4.


texts:
I have a large polygone file with landuse information (6 classes). Its vectorised from a raster classification so the polygons are all angular. Any of my tries to generalize with v.generalize (sneaks, chaiken, douglas...) ends after minutes of working in a programm crash (GRASS and Qgis using sextante).

How can I smooth the edges of my polygons? Thanks!


texts:
I've try to convert SHP files to geoJSON. When batch process is started occure problem. Here is copy flrom log:

Processing algorithm 1/233…
Algorithm Convert format starting…
Input parameters:
{'INPUT': 'sp gid_1', 'OPTIONS': '', 'OUTPUT': 'D:/test/1sp gid_1.geojson'}

GDAL command:
ogr2ogr -f "GeoJSON" "\"D:/test/1sp gid_1.geojson\"" C:/Users/vimars/AppData/Local/Temp/processing_01832e3931f144a18ec42d0c69d8bdb5/fb63b857189c43f9b5180a5c6b4b80ed/INPUT.shp gid_1
GDAL command output:
ERROR 4: Failed to create GeoJSON datasource: "D:/test/1sp gid_1.geojson":

ERROR 1: GeoJSON driver failed to create "D:/test/1sp gid_1.geojson"

Algorithm Convert format correctly executed…
Execution completed in 0.36 seconds
Results:
{'OUTPUT': 'D:/test/1sp gid_1.geojson'}

Loading resulting layers

I've triple chech that d:\test have write permission but still no luck :( No file was created in destination folder.

Additional info: single task work with no problem:

Processing algorithm…
AlgorithmDialog
Input parameters:
{ 'INPUT' : 'D:\sp|layername=gid_1', 'OPTIONS' : '', 'OUTPUT' : 'D:/test/1.geojson' }

GDAL command:
ogr2ogr -f "GeoJSON" "D:/test/1.geojson" C:/Users/vimars/AppData/Local/Temp/processing_01832e3931f144a18ec42d0c69d8bdb5/cd83614a8607490c9e6eb526af0c309c/INPUT.shp INPUT
GDAL command output:
Execution completed in 0.85 seconds
Results:
{'OUTPUT': <QgsProcessingOutputLayerDefinition {'sink':D:/test/1.geojson, 'createOptions': {'fileEncoding': 'System'}}>}

Loading resulting layers
Algorithm 'Convert format' finished


texts: How to show just a few objects in the attribute table?

texts:
I am pretty lost with the attribute table, so I couldn’t figure out the already existing answers about it (using QGIS 2.12.3.).
I have a shapefile with all the protected areas in one region. Since I don’t need all of those protected areas, I would like to know how I can hide all the other objects?
I have to search for the objects through their name and really don’t know how to use the dialog with all its functions.


texts: Python GDAL RasterizeLayer option parameter

texts:
I'm writing some code that's supposed to take multiple multipolygon layers and merge them into 1 raster. I can load in the datasets just fine and render them to an image (although I use this just for debugging), here's what I'm doing right now:

noiselayers["noise"] = [[1, 0.1], [2, 0.8]]

    #create a source in memory
    ds = ogr.GetDriverByName("Memory").CreateDataSource("wrk")
    source = ds.CreateLayer("", None, ogr.wkbMultiPolygon )
    source.CreateField(ogr.FieldDefn("VAL", ogr.OFSTFloat32))

    for layerdata in noiselayers[selected]:
        data = requests.get("http://localhost/api/layer/export/%s" % str(layerdata[0]))
        geometry = json.loads(data.text)

        for g in geometry:
            geom = ogr.CreateGeometryFromWkt(g['the_geom'])

            structure = ogr.Feature(source.GetLayerDefn())
            structure.SetGeometry(geom)
            structure.SetField("VAL", layerdata[1])

            #save geometry to layer
            source.CreateFeature(structure)

    pixelsize = 1000
    #novalue = -9999

    #get the layer bounds dynamically from the layer
    #x_min, x_max, y_min, y_max = layer.GetExtent()

    #extracted values from ns_play_area_final_EPSG_3035
    x_min = 2777236.78057
    x_max = 5065039.36824
    y_min = 2815620.0976
    y_max = 4485652.62185

    x_res = int((x_max - x_min) / pixelsize)
    y_res = int((y_max - y_min) / pixelsize)

    target = gdal.GetDriverByName("MEM").Create("", x_res, y_res, 1, gdal.GDT_Byte)
    target.SetGeoTransform((x_min, pixelsize, 0, y_max, 0, -pixelsize))

    gdal.RasterizeLayer(target, [1], source)

This grabs a bunch of data from our API and rasterizes it just fine. However, I need to merge multiple layers with different weights (in the noiselayers variable, the 2nd value of each list is the weight). When it's a higher value it should increase the value more, for now just setting the value of the square to that value is enough. I've tried using the options (as seen here: https://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py), but they never seem to do anything, here's some of the things I tried to use the options:

gdal.RasterizeLayer(target, [1], source, None, None, [1000], ["ATTRIBUTE=VAL"])

gdal.RasterizeLayer(target, [1], source, options = ["ATTRIBUTE=VAL"])

any combination of the 2, the options never seem to do anything other than break the rasterizing. I've also tried inversing (using INVERSE=TRUE, or any combination of capitalizing/not capitalizing that) the raster, but that also never seemed to work.

Does anyone have a working example of how to use attributes and/or other hints on how to solve my issue?


texts:
As title says , I've tried to look for it and didn't find a way to do it.
Did I miss it or it wasn't implemented yet ?

I consider switching to OL 3 because of the export feature (I want to add easy print of application).


texts: GRASS GIS r.in.gdal, import multiple band ncdf

texts:
I am trying to import a ncdf file with 20000 bands (time steps of precipitation) to grass gis using r.in.gdal. My problem is that i dont need all the bands, only a portion of it, e.g : band 1 to band 200, then band 200 to band 400, etc. The only way i can do that is manually put the numbers of the band in module config '[band=integer[,integer,...]] ', e.g band = 1,2,3,4,5... 200.

Is there a way that a i can put for example like in python 'range(1,200)' and get the 200 first bands of the file?.


texts:
I am new to QGIS. I just happened to browse to a tutorial on georeferencing in QGIS and I found it easy and interesting.

I also calibrated the map and saved the points as per WSG84 standard. Now after I closed the software and when I restarted I am unable to see the coordinates when I load the GeoTIFF image.

Why am I not able to load the coordinates?


texts:
For web-mapping with OpenLayers 3x and JavaScript (front-end), is it possible to exactly clip ("made to look") images/rasters (Landsat-8 @ 15 & 30 m, Sentinel-2 @ 10 m) to the edges of irregular polys? Clipped perfectly by CSS or alpha masking or other approaches?
I guess I can resample to 5 m or something, however I have other functions where I allow the user to click to query pixel values and I don't want to mislead the user to what they click to against what they see..


texts: Editing live MSSQL Tables and Updating Fields

texts:
I am looking at creating a trigger on a live MS SQL spatial table that will update a field when a user either:

Creates a new feature
Updates the geometry of a feature.

Unfortunately I have only been able to work out updating a field when a new feature is created, or updating a field for every record when the geometry of a record is changed.

Ideal world - attribute 'Updated' populates current date for the particular feature that has been created or modified.

This is what I have so far that works only when creating a new feature:

SET ANSI_NULLS ON
GO
SET QUOTED_IDENTIFIER ON
GO

CREATE TRIGGER *[triggername]*
ON  *[table]*
AFTER INSERT,UPDATE
AS
BEGIN
DECLARE @ID VARCHAR(8)
SELECT @ID = [ID] FROM *[table]*
UPDATE *[table]*
   SET Updated = GETDATE()
   WHERE [ID] = @ID
END
GO


texts: Add feature to an OpenLayers VectorTileSource

texts:
Is there any way I could add a Feature to an OpenLayers VectorTileSource?
I know that the VectorSource has the method addFeatures() but not the VectorTileSource since it's of type MVT.
What I need to do is to get the Features existing in a VectorTileSource and add a dummy feature to the collection for a small test. If I could find a way of doing it would be great.


texts:
I am using QGIS and I have a raster from which I would like to extract all the pixel values and if possible their corresponding coordinates.

I am able to use the identify feature option to point to a location on the raster to get the individual raster value; but, I would like to do this for all points without having to click all of them.


texts: QGIS Freehand editing plug-in: blocked attribute table

texts:
I can't save object after drawing by Freehand editing plug-in. I suppose it will be connected with creating some (special) column in the table.

Could anyone suggest a solution?


texts: im trying to install the toolbox Swath profile, but the program send me an error "NameError: name 'Line_shapefile' is not defined". How can i install it?


texts: Can gdal2tiles create tiles at zoom levels > 10?

texts:
I am trying to create tiles from tif images using the gdal2tiles command (at the OSGeo4W command line). I've been asked if it would be possible to create tiles at zoom level 12. Using the -z switch like so ...

gdal2tiles --s_srs EPSG:32717 -z 12 input_tif.tif output_folder

This creates the tiles but they are all black. I suspect this cannot be done but I'd like to confirm that I'm not just having an unrelated problem. Tiles created at 1-10 appear to be ok with most of my input files.


texts: Entering UTM coordinates with decimals using ArcGIS Pro

texts:
I have UTM coordinates containing decimals (e.g. 884305.145mN). I need to create a polygon using the Absolute X,Y,Z feature in ArcGIS Pro.
But ArcGIS Pro doesn't seem to allow decimals in UTM degrees.
How do I make ArcGIS Pro accept UTM coordinates that contain degrees?


texts: Running ArcMap Addin Button with VBA from MS Access form?

texts:
I am using ArcGIS Desktop 10.0.

I need a way to execute an ArcGis-Addin (vb.net) from a MS-Access database form (VBA).

The problem is that there is no ESRI.ArcGIS.Desktop.Addin reference in VBA and the "pItem.Execute" method with the "pCmdBars.Find(pUID)" does not work with my Addin-Button.

I want to open the ArcMap.mxd file and perform a selection that shows the spatial data of my MS-Access database query.

The rest is not a problem  (maybe!).


texts: Qgis simple line vertex geometry and Z values

texts:
QGIS 2.8.4 - utility networks
postgres 9.3/postgis 2 database

The linear network is made up of simple lines joining points that have different z values. For each feature in the line layer's table I would like there to be a column for the start point and a column for the end point populated with the z values from the point layers.

I would like to update a column in a the line attribute table where the first vertex of a simple line is updated with a numerical value from another point layer to which it is snapped.
Then do the same for the second vertex and get the value from the point later to which the second line point is snapped.

From there a simple calcule to get the slope.for each line in the network.

Would ST_startpoint and ST_endpoint be the right way to go?


texts: Why doesn't QGIS dissolve work properly for some datasets?

texts:
I am used to work with ArcGIS and now I started to explore QGIS more. I have problem with dissolve function. I have administrative borders and I would like to dissolve them into higher administrative areas. When I process dissolve according to common attribute NUTS3 - I get strange results - see pic.

It works fine in ArcMap 10.2.2. can you tell me, what do I do wrong in QGIS?


texts:
QGIS Version: 2.14.3
GRASS Version: 7.2

Both installed via Homebrew.

Tried setting the GRASS7 path in "Processing->Options" but QGIS does not accept it.

Related to the following questions:

QGIS 2.14.1 LAStools install error: "wrong value for parameter 'Msys folder'..."

QGIS GRASS Missing Dependancy

Grass7 processing toolbox QGIS

Update:

I modified grassPath() in "Grass7Utils.py" and changed the grass path folder setting to:

folder = '/usr/local/Cellar/grass-72/HEAD-/grass-7.2.svn/'
return folder

Checking the processing options in QGIS, the path change was accepted. However, GRASS7 algorithms still come up with the error message that GRASS7 is not configured correctly. Which GRASS path does QGIS need?


texts: Missing / misnamed fields in PostgreSQL table join in QGIS

texts:
First, I'm brand new to PostgreSQL and PostGIS, so my question is probably asked wrong, and it may be a duplicate – but my search attempts didn't find anything. Also, I don't expect this is an issue with the data_ ¯_(ツ)_/¯ _but, for reproduction reasons, the tables are the 2020 and 2019 CSVs from Boston's property assessment.
Problem: When I join two tables with matching field names in the QGIS DB Manager, only one set of fields appear; when I rename the fields, the wrong fields are renamed. When I run the same commands in a pgAdmin window, these problems don't happen.
What I did & what happened
I tried joining two tables in the QGIS DB Manager
SELECT
    a.pid, a.owner, b.pid, b.owner
FROM
    properties_20 AS a
        INNER JOIN properties_19 AS b ON a.pid = b.pid
LIMIT 3;

Resulting view:
pid         owner
0502309029  LADAKH REALTY LLC MASS LLC
0502309031  ROGERS ALEX A
0502309032  CLAIRE A WALTON 2012

The resulting view only had two fields: pid and owner. I expected four: a pid for each table and an owner for each table. As I'm new to all of this, though, I thought it might be a matter of matching field names, so I renamed the owner fields with AS.
SELECT
    a.pid, a.owner AS owner_20, b.pid, b.owner AS owner_19
FROM
    properties_20 AS a
        INNER JOIN properties_19 AS b ON a.pid = b.pid
LIMIT 3;

The resulting view added a third column named owner_19, but the column's contents were the table's pid values:
pid         owner_20                     owner_19
0502309029  LADAKH REALTY LLC MASS LLC   0502309029
0502309031  ROGERS ALEX A                0502309031
0502309032  CLAIRE A WALTON 2012         0502309032

I then tested the same commands outside of QGIS, using pgAdmin 4:
First query (no renaming)
pid             owner                           pid             owner
"0502309029"    "LADAKH REALTY LLC MASS LLC"    "0502309029"    "LADAKH REALTY LLC"
"0502309031"    "ROGERS ALEX A"                 "0502309031"    "ROGERS ALEX A"
"0502309032"    "CLAIRE A WALTON 2012"          "0502309032"    "CLAIRE A WALTON 2012"

Second query (renaming)
pid             owner_20                        pid             owner_19
"0502309029"    "LADAKH REALTY LLC MASS LLC"    "0502309029"    "LADAKH REALTY LLC"
"0502309031"    "ROGERS ALEX A"                 "0502309031"    "ROGERS ALEX A"
"0502309032"    "CLAIRE A WALTON 2012"          "0502309032"    "CLAIRE A WALTON 2012"

Primary question: Why is this happening? Is there a way I can fix this so it works in QGIS? Is my syntax somehow lazy?
Secondary question: If I should just be working in the pgAdmin window, how do load the results as a new layer in QGIS?
I considered posting this on the StackOverflow site, but because the problem is specific to QGIS I asked here.


texts: ESRI ArcGIS error: Zero length polylines not allowed

texts:
While writing a simple ICommand in ArcMap 10.2.1 using C# that lets you select two different point features on the map and it creates a feature in a polyline feature class with those two points as its end points (connect the dots!) using the simple logic:
grab the two selected points
add the two points to an IPointCollection
call IFeatureClass.CreateFeature() on the polyline feature class
set the new feature's shape to the point collection
call IFeature.Store() to save it.
Then we got this error: "Zero-length polylines not allowed..."
Even adding a few debug lines in to check the new feature's shape correctly has the FromPoint and ToPoint coordinates AND the Shape.Length correctly shows 204.34566 it still gives this error.


texts:
I'd like to develop an app so that user introduce the layers in comboboxes ( layers are slope,aspect,elevation etc) then user use these layers for processing in another form(processing form) and  dose not need to select the items in processing form . i don't know how can i write the code for saving the layers in comboboxes .

Can you recommend a solution for this app?


texts: Why is Orthorectify button not active in Image Analysis Windows in ArcGIS10?

texts: I wanted to Ortho rectify GeoEye HR imagery using existing RPC and IFSAR DTM in Image Analysis. According to ArcGIS Resource center http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/using_the_Orthorectify_button/009t000000mv000000/ the button is automatically active if the imagery has RPC file. In my case the ortho button was not active although my Imagery has RPC file. Please advice thank you.


texts: Does not display Map in Rstudio using Python Package - GEEMAP

texts:
I am trying to use the GEEMAP package of Python in RStudio-1.4.1106.
After importing the packages and running the code shown below MAP is not displayed rather it displays the MAP attributes in the console:
Map(center = [40, -100],
  controls = (
    WidgetControl(
      options = ['position'],
      widget = HBox(
        children = (
          ToggleButton(
            value = False,
            icon = 'globe',
            layout = Layout(
              height = '28px',
              padding = '0px 0px 0px 4px',
              width = '28px'),
            tooltip = 'Search location/data'), ))),
    ZoomControl(options = ['position', 'zoom_in_text',
      'zoom_in_title', 'zoom_out_text', 'zoom_out_title'
    ]),
    ScaleControl(options = ['imperial', 'max_width', 'metric',
      'position', 'update_when_idle'
    ], position = 'bottomleft'),
    FullScreenControl(options = ['position']),
    MeasureControl(
      active_color = 'orange',
      options = ['active_color', 'capture_z_index', 'completed_color', 'popup_options', 'position',
        'primary_area_unit', 'primary_length_unit', 'secondary_area_unit', 'secondary_length_unit'
      ],
      popup_options = {
        'className': 'leaflet-measure-resultpopup',
        'autoPanPadding': [10, 10]
      },
      position = 'bottomleft', primary_length_unit = 'kilometers'),
    AttributionControl(
      options = ['position', 'prefix'], position = 'bottomright'),
    DrawControl(
      circle = {
        'shapeOptions': {
          'color': '#3388ff'
        }
      },
      marker = {
        'shapeOptions': {
          'color': '#3388ff'
        }
      },
      options = ['position'],
      polygon = {
        'shapeOptions': {}
      },
      polyline = {
        'shapeOptions': {}
      },
      rectangle = {
        'shapeOptions': {
          'color': '#3388ff'
        }
      }),
    WidgetControl(
      options = ['position'],
      position = 'topright',
      widget = VBox(
        children = (
          ToggleButton(
            value = False,
            icon = 'wrench',
            layout = Layout(
              height = '28px',
              padding = '0px 0px 0px 4px',
              width = '28px'),
            tooltip = 'Toolbar'), )))),
  crs = {
    'name': 'EPSG3857',
    'custom': False
  },
  default_style = MapStyle(),
  dragging_style = MapStyle(cursor = 'move'),
  layers = (
    TileLayer(
      base = True,
      max_zoom = 19,
      min_zoom = 1,
      options = ['attribution', 'detect_retina', 'max_native_zoom', 'max_zoom', 'min_native_zoom', 'min_zoom', 'no_wrap', 'tile_size', 'tms']),
    TileLayer(
      attribution = 'Google',
      name = 'Google Maps',
      options = ['attribution', 'detect_retina', 'max_native_zoom', 'max_zoom',
        'min_native_zoom', 'min_zoom', 'no_wrap', 'tile_size', 'tms'
      ],
      url = 'https://mt1.google.com/vt/lyrs=m&x={x}&y={y}&z={z}')),
  layout = Layout(height = '600px'),
  options = ['bounce_at_zoom_limits', 'box_zoom', 'center', 'close_popup_on_click', 'double_click_zoom', 'dragging', 'fullscreen', 'inertia', 'inertia_deceleration', 'inertia_max_speed', 'interpolation', 'keyboard',
    'keyboard_pan_offset', 'keyboard_zoom_offset', 'max_zoom', 'min_zoom', 'scroll_wheel_zoom', 'tap', 'tap_tolerance', 'touch_zoom', 'world_copy_jump', 'zoom', 'zoom_animation_threshold', 'zoom_delta', 'zoom_snap', 'zoom_start'
  ],
  scroll_wheel_zoom = True,
  style = MapStyle(), zoom = 4.0)

import numpy as np
import geemap
import ee
Map = geemap.Map(center=(40, -100), zoom=4)
Map

How can I display the map?


texts:
I am a new user of ArcGIS online. I need to build a map that allows to see the name of the inspector in charge of the area.
So far I succeeded in creation a map that shows the boundaries of the neighborhoods and a pop-up (very basic) that shows the information required.
I would like to know if there is a way to show the same pop-up (for the corresponding area) when I use the address locator.
Example: If I type the address of my house I would like to see (right away) a pop-up showing the name of the inspector that I should contact.


texts: Chaining outputs in ArcPy standalone script

texts:
I am trying to create a standalone ArcPy script.
I started by creating this in ModelBuilder and then I exported it as a Python script. This works perfectly in ModelBuilder but produces only an empty Shapefile in the standalone script!
The script is pretty simple:

Create a layer out of 2 enterprise geodatabase as input to 'Add_Join' tool
Join the two layers
Create a view out of a standalone table in the same geodatabase
join the table view to the result of step 2
use 'FeatureClassToFeatureClass_conversion' to save as shapefile

For some reason I get only an empty shapefile. The final 'FeatureClassToFeatureClass_conversion' tool does not seem to wait for the joins to complete.
Is a parameter dependency an option here for standalone ArcPy tools?
Primary_FeatClass = "C:\\Temp\\PROD_GIS_XS01_NIS.sde\\SchemaName.Electric\\SchemaName.Primary_FeatClass"
FeatClass2 = "C:\\Temp\\PROD_GIS_XS01_NIS.sde\\SchemaName.Electric\\SchemaName.FeatClass2"
Table2 = = "C:\\Temp\\PROD_GIS_XS01_NIS.sde\\SchemaName.Table2"

Field_Map = "INSTALDT \"Installed Date\" true true false 36 Date 0 0 ,First,....<more text>...." << truncated

# ----------------Join 2 feature layers---------------
#Feature layer from GDB feature class
arcpy.MakeFeatureLayer_management(Primary_FeatClass, Primary_FeatClass_Layer, "", "", "OBJECTID OBJECTID VISIBLE....<more text>....") << truncated

#Feature layer from GDB feature class
arcpy.MakeFeatureLayer_management(FeatClass2, FeatClass2_layer, "", "", "OBJECTID OBJECTID VISIBLE NONE;....<more text>....") << truncated

#Add Join them together to create the first_join_result
arcpy.AddJoin_management(Primary_FeatClass_Layer, "RELSUPPORTSTRUCTUREOBJECTID", FeatClass2_layer, "OBJECTID", "KEEP_COMMON")

# ---------------- Join another Table to the above-----------------------
#Make Table View
arcpy.MakeTableView_management(Table2, Table2_View, "", "", "OBJECTID OBJECTID VISIBLE NONE;....<more text>....") << truncated

#Join table view to the first_join_result
arcpy.AddJoin_management(Primary_FeatClass_Layer, "SchemaName.Primary_FeatClass.CISSERVICEACCOUNTNUMBER", Table2_View, "SERVICEACCOUNTASSOCIATIONCD", "KEEP_COMMON")

#--------------- export to Shapefile ----------------------

# Process: Feature Class to Feature Class
arcpy.FeatureClassToFeatureClass_conversion(Primary_FeatClass_Layer, TEMP, "Primary_FeatClass.shp", "", Field_Map, "")


texts: Populating Dropdown based on Geoserver service

texts:
I am new to web programming here. I have a table which has multiple columns including: Restaurant Name,Restaurant Id, Lat, Lon, City, District, Country.

My requirement is to build a web client using OpenLayers and JavaScript and allow users to search a Restaurant. I am using GeoServer to create a layer from this table.

My problem is : I have to add dropdowns i.e. one for country, another for District /City based on the country selected. This city/district/country mapping has to come from the table which is in Teradata. Is there any way I can get the data to populate the dropdown from GeoServer services Or do I need to create a new service of my own.
If I create my own service, can it be hosted from GeoServer.

Any inputs?


texts: Can I get time & distance from OSM2PO without loading data into PostgreSQL?

texts: Can I get time and distance from OSM2PO routing without loading the data into PostgreSQL? I've seen the answer regarding matching IDs with the DB, but I need to do this using only the OSM2PO service.


texts: Looking for a raster calculator C# class

texts:
I am developing an application using C# and dotspatial. I need to performs raster calculation a lot..I need a library which helps me to do it..Do you know any library for this purpose?
The library must consider boundary of input rasters and also input bands and pixel sizes..

I only need some operations such as  +,-,* and /


texts: Create Layer Package fails to execute -- error code 00045

texts:
I'm experiencing an error when I attempt to save a feature class as a Layer Package (.lpk) by right-clicking that feature class in my Table of Contents and selecting "Save As Layer File".

The Layer Package form pops up. After entering the required bits (namely summary and tags) and clicking "Analyze", I receive an error (code 00045) stating:

  A layer description is required for packaging

First, this field is not shown as required within the Layer Package form. Second and more importantly, even after entering a description, the error continues and prevents me from creating a Layer Package.

No amount of fiddling with the description will allow the layer package to be created. I can find no way around it. Esri's help page for code 00045 offers no assistance in remediating this particular problem, nor do I see that others are experiencing or have experienced it.

Can you suggest a remedy?


texts: How to perform multiple string transformations using one transformer in FME Workbench?

texts:
I am converting a .DWG file to .GDB. All fields in DWG file are of Char type. The target .GDB needs to have some attributes in Float data type. For now, I am using one StringFormatter transformer for conversion of one attribute. That is cluttering the workspace.

Can we perform all Char to Float transformations using a single transformer?

Thanks


texts: Polygonize area surrounded by buffer in QGIS

texts:
I have an automated system in QGIS that separates loops from GPX traces and works out their centroid using 'Polygonize' >> 'Keep N biggest parts' >> 'Centroids'.
However, quite often it includes part of the 'loop approach' in the polygon:

This distorts the centroid calculation.
I think there might be a way to use a line buffer to separate the loop from the GPX trace. I don't care about preserving the GPX data but I want to avoid 'dangles'.
I can use a buffer to make a polygon around the GPX trace:

But I cannot work out how to make the area inside the 'buffer polygon' (shown in yellow) into a separate polygon that I can use to calculate the centroid.

It is essential that the process can be automated. I have hundreds of files that need to be processed.
Negative buffers don't help:

They change the buffer but don't make a polygon from the area surrounded by the buffered line.
But maybe I am missing something.
here is the GPX trace for the simple square so you can test it:
<?xml version="1.0" encoding="UTF-8"?>
<gpx xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.topografix.com/GPX/1/1" xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd" version="1.1" creator="Open GPX Tracker for iOS">

        <trk>
        <trkseg>
            <trkpt lat="100.000" lon="100.000">

            </trkpt>

            <trkpt lat="100.000" lon="200.000">

            </trkpt>

            <trkpt lat="200.000" lon="200.000">

            </trkpt>

            <trkpt lat="200.000" lon="100.000">

            </trkpt>

            <trkpt lat="100.000" lon="100.000">

            </trkpt>

        </trkseg>
    </trk>
</gpx>


texts: why use geoserver layername workspace prefix?

texts:
I am using openlayers to consume geoserver. I have created a geoserver WMS layer and preview it.

var tiled = new OpenLayers.Layer.WMS(
                    "GIS:pipelines - Tiled", "http://10.1.1.9:8080/geoserver/GIS/wms",
                    {
                        LAYERS: 'GIS:pipelines',
                        STYLES: ''
                    },
                    {
                        isBaseLayer: false
                    }
                );

I have two questions:

http://10.1.1.9:8080/geoserver/GIS/wms Why URI is containing GIS workspace name.
LAYERS: 'GIS:pipelines' is using GIS workspage prefix.

If I delete GIS from URI and LayerName, may be error?


texts: Getting XY coordinates and field values for each row in ArcPy SearchCursor?

texts:
I am trying to enumerate each row of a points shapefile and get a certain field's value (elevation) and the xy coordinates.

But I see I can only use:

arcpy.SearchCursor to get the field; OR
arcpy.da.SearchCursor(fc, "SHAPE@XY") for the xy coordinates

Can I do both at once?


texts: How to programmatically load an annotationlayer to ArcMap

texts:
I have a geodatabase with an annotation feature class in it containing 3 subtypes. Now I load the featuredataset to the map with this code:

IFeatureClass featClass = (workspace as IFeatureWorkspace).OpenFeatureClass(sNameOfFeatclass);
IFeatureLayer fLyr = new FeatureLayer();
fLyr.FeatureClass = featClass;

In ArcMap the opened Layer is displayed but symbology is like polygon classified with the subtypes. In debug mode I see that the featureclass is of type esriFTAnnotation and the subtypes are present but the layer does not support the interface IAnnotationLayer.

In resources I only have found how to create an AnnotationFeature Class, but in my case the featureclass exists. How can I load it to map so that in layer properties the tab Annotation is displayed and not Labeling?


texts: Getting img.id is not a function error when calculating VCI

texts:
I try to calculate VCI, when I run this code I got an error

  img.id is not a function

How to fix this problem?

//Image listing
var modis = ee.ImageCollection("MODIS/006/MOD13Q1").filterBounds(table).select('NDVI');

var modis_list = [];
for(var i = 6; i <=10; i++) {
  for(var j = 2000; j <=2019; j++){
    modis_list.push(modis.filter(ee.Filter.calendarRange(i, i, 'month'))
                          .filter(ee.Filter.calendarRange(j, j, 'year')));
  }
}
print(modis_list);

//VCI calculation
var vci = modis_list.map(function(img){
 var id = img.id();
 var min =  img.reduceRegion(ee.Reducer.min(), geometry,250).get('NDVI');
 var max = img.reduceRegion(ee.Reducer.max(), geometry,250).get('NDVI');
 return img.expression(
   "(NDVI-min)/(max-min)",{
     "NDVI" : img,
     "max" : ee.Number(max),
     "min" : ee.Number(min)
   }).copyProperties(img,['system:time_start','system:time_end']);
});
print(vci)

I have no prior knowledge about the Google earth engine, I just try to learn.


texts: Multiline VBScript Label Expressions in Arcpy

texts:
I am able to set simple label expressions, but when I try to set multiline expressions (which all VBScript expressions are) it does not get copied into the expression box properly.  I have tried using Python's triple quote technique for multilines but that only seperates each line with a tab in the expression window.  Anyone know of a way to do this?  If so, does the label class need to be set to 'advanced' first, or is it possible to automate this aspect as well?

Alternatively, loading in precreated expression files would work just as well if not better if this is possible to script.

This is with the use of the Maplex extension.

EDIT

Basically this code below needs to be pasted into the expression box with Arcpy.  One line works just fine, but because there are multiple lines, I can't get it to paste it properly.

Function FindLabel ( [SUPP_NO], [SUPP_SUFF])
IF [SUPP_SUFF] = "000" THEN
FindLabel = "#" & [SUPP_NO] & VBNewLine & [Fuse_Count]"


texts: Programmatically selecting features by polygon

texts: Expanding upon this post how could you use QGIS to programmatically draw a selection polygon and select features within the polygon?  I have been using QgsMapToolEmitPoint(self.canvas), to find and select features by point but I am struggling to find a similar polygon solution.


texts: ArcGIS how to geolocate/georeference raster data arrays?

texts:
How does one geolocate (or georeference) a raster of swath grided data with separate lat/lon arrays and multiple data bands using ArcMap 10.3?

A sample data file (*.h5) can be downloaded from here

A simplified data structure is:

Subdataset 1: lats

Subdataset 2: lons

Subdataset 3 to end : data bands


texts: How to merge (append) two tables in ArcGIS 9.3.1?

texts: How do I merge 2 tables (they both have the same field names and columns) in ArcGIS 9.3.1 ?


texts: Removing specific Circle Markers in Leaflet Shiny

texts:
I have a layer of CircleMarkers and I am trying to remove only the markers that have a certain layerId. The id's for these circle markers are in a dataframe. In reality I have millions of id's whose subset to visualize depends on a different dataset which can change over time. So I need a dynamic way to delete specific layerId's.

Below is a simple example:
Suppose I have a dataframe with 3 rows with id's 1, 2 and 3. I tried to make a checkboxInput with the options to delete id's 1 and 2 or 3.

Below the inputs will trigger an ObserveEvent that use the removeMarker function. However, nothing happens. I have tried a million ways to enter the id's into the removeMarker and I have also tried several of the other ways to deletion. Either nothing happens or all disappear. I need a way to delete specific markers. Note that this in no way reflects my real life example: I have millions of points and deletion of points is not through a direct checkbox input, but based on a different dataset that can change over time.

 ui <- shinyUI(fluidPage(
sidebarLayout(
    sidebarPanel(
        checkboxInput("delete1", "Delete ID=1 and 2",value=FALSE),
    checkboxInput("delete3", "Delete ID=3",value=FALSE)
    ),
    mainPanel(
        leafletOutput("map")
    )
)
))

df <- data.frame(id=c(1,2,3),lng = rnorm(3, -106.1039361, 0.5) ,
              lat = rnorm(3, 50.543981, 0.5))

server <- shinyServer(function(input, output, session) {

output$map <- renderLeaflet(
    leaflet() %>%
addTiles() %>% addCircleMarkers(layerId=df$id,df$lng,df$lat, group='marker', radius=2, fill = TRUE,color='red')

    )

observeEvent(input$delete1, {
    proxy <- leafletProxy('map')
    if (input$delete1){ proxy %>% removeMarker(df[1:2,1])
 }
 })

observeEvent(input$delete3, {
    proxy <- leafletProxy('map')
    if (input$delete3){ proxy %>% removeMarker(3)}
   })
})

 shinyApp(ui, server)


texts: FME: how to replace a character at a specific index?

texts:
Suppose I've a string foo,
How can I change the 3rd character (identified by its index) ?
The Ruby equivalent would be
# Define the variable string and assign the value 'bar'
string    = 'foo'     #=> 'foo'

# Change the string's character at index 2 to 'Z'
string[2] = 'Z'

# Output the string variable to see changes
string                #=> 'foZ'

I've tried using SubString(@Value(my_string), 2, 1) with replacement text set to 'Z', but it outputs 'fZZ' (replacing all occurences of the value of the character at the specified index).
I could archive it by using regular expressions, though it'd be an overkill


texts: Display ESRI MapService (WMTS) with EPSG:28992 coordinate system in OpenLayers web map with OpenStreetMap base layer

texts: ESRI MapService https://tiles.arcgis.com/tiles/nSZVuSZjHpEZZbRo/arcgis/rest/services/Historische_tijdreis_1870/MapServer have EPSG:28992 coordinate system. In QGIS I can add this layer via WMTS https://tiles.arcgis.com/tiles/nSZVuSZjHpEZZbRo/arcgis/rest/services/Historische_tijdreis_1870/MapServer/WMTS (layer: Historische_tijdreis_1870, format: image/jpgpng, tileset: default028mm, crs: EPSG:28992). How I can display this layer in OpenLayers (6) web map?


texts: Map opening in Indian Ocean on QGIS Web Client?

texts: After installing web part of qgis i created a project with layers. All layers are working, but i always need to scroll the map to reach these layers. How can I open the map exactly there i want? I tried to set epsg in file GlobalOptions.js, created file in the path site/libs/proj4js/defs/ with the correct arguments, added epsg parameter to openlayers.js, but still it opened indian ocean. Also in the project under the crs properties i added right crs.


texts: How can one create a binary image from a FeatureCollection in Google Earth Engine?


texts: Can I force an attribute in my WFS Layer to be unique, thus, don't allow mutiple feature with the same value for the attribute?


texts: When I save a .shp file as .kml and keep the DescriptionField as "Description", the Description field in the attribute table is blank, but when the DescriptionField is set to "Descriptio" (no "n" character on the end), the Description field is properly completed.


texts:
i have 2 layers. layer1 and layer2. Both are polygons. I want to get both layers Shape_Area. I know I can get there shapeareas by using a searchcursor and extracting them, which I have done already.

I want to know if there is a quicker, shorter way to get a layers shapearea than using two separate search cursors?

My long way:

areas = {}
with arcpy.da.SearchCursor(layer1, "Shape_Area") as cursor:
    for x in cursor:
        areas['layer1'] = x
with arcpy.da.SearchCursor(layer2, "Shape_Area") as cursor2:
    for y in cursor2:
        areas['layer2'] = y

I know there is arcpy.Polygon.getArea() but I am unsure how to use it...


texts: How to remove overlaps and keep highest priority polygon

texts:
I have a multipolygon table with overlapping and non overlapping multipolygons. Where they overlap I need to keep the one with the highest (=1) priority and erase away from the lower priority ones to make the table non-overlapping:

create table table1 (
prio    int,
geom    geometry(MULTIPOLYGON,3006));

insert into table2 (prio, geom) values
    (1, ST_PolygonFromText('MULTIPOLYGON(((666263 7266050,666869 7259810,671789 7265158,666263 7266050)))')),
    (2, ST_PolygonFromText('MULTIPOLYGON(((670720 7261236,668652 7264445,672966 7264124,670720 7261236)))')),
    (3, ST_PolygonFromText('MULTIPOLYGON(((669401 7267012,670328 7263233,674677 7266157,669401 7267012)))')),
    (4, ST_PolygonFromText('MULTIPOLYGON(((666263 7258563,666905 7254748,669615 7257315,666263 7258563)))'));

I can maybe solve this with something like this but its going to require alot of code since I have ~40 prios.
select id, prio, geom from table1 where prio=1
union
select p1.id, p1.prio, ST_Difference(p2.geom, p1.geom) geom from
(select * from table1 where prio=1) p1
join
(select * from table1 where prio=2) p2
on st_intersects(p1.geom, p2.geom) and p1.id <> p2.id

--union
--select id, prio, geom from table1 where prio=2
-- ...

Isnt there some easier way?
(The database should have been created with a constraint(?) that doesnt allow overlap, but it wasnt... Now I have ~45k polygons that partially overlap)


texts: Configuring Mapproxy layer with non-ansi characters

texts:
I have a working Mapproxy 1.6.0 installation serving several WMS layers.
I am now tasked with adding an extra WMS-layer which contain the Danish character Ø.

I have defined the layer in Mapproxy.yaml like this:

layers:
  - name: Other-objects
    title: Other objects
    sources: ['wms_trans:Øvrige_objekter']

When starting Mapproxy it quits immediately with the following information in the log:

2015-03-11 09:42:30,322 - INFO - mapproxy.config - reading: c:\mapproxy\mapproxy.yaml
2015-03-11 09:42:30,338 - WARNING - mapproxy.config - u'wms_trans:\xd8vrige_objekter' in layers[3].sources[0] not of type str
2015-03-11 09:42:30,338 - CRITICAL - mapproxy.config - invalid configuration

The Ø have been translated into \xd8.
I have attempted to save mapproxy.yaml with ansi, utf8 and utf8-signature encodings, this does not seem to have an effect.

Since this is an external wms source, it is not possible to rename it.
How can I make this work?


texts: Unable to activate GDAL Georeferencer in QGIS

texts:
I have installed QGIS 2.14 on my MacBook Air, but am unable to activate the "Georeferencer GDAL" plugin.  When I open the "Manage and Install Plugins" dialog, the Georeferencer GDAL plugin is listed under Installed, but it is not activated.  I tried clicking on the checkbox for the plugin, but that has no effect.  Same for all of the other plugins listed as installed but not activated.  (Also, the uninstall and reinstall buttons are disabled.)

How can I activate georeferencing?


texts: Efficiently identifying overlapping polygons with Python

texts:
I'm parsing the NOAA severe weather alert feed and using it as a datasource for a weather alerting system.

In the feed, I occasionally find a polygon entry for an alert. Unfortunately, counties that fall under these polygons are not always identified in the FIPS/UGC codes section of the alert.

In cases where I come across a polygon entry, I need to interate over all my county polygons to make sure all counties potentially affected by weather are included in my output.

Here's an example of a polygon entry from today:

<cap:polygon>41.38,-97.93 41.40,-97.81 41.45,-97.69 41.40,-97.69 41.36,-97.92 41.38,-97.93</cap:polygon>

With 3200+ US County polygons that I have stored, I need to find the most efficient way of iterating over them and comparing the county polygon with the alert polygon to see if they overlap. (I'm not sure if I'm looking for overlap, intersection, or union?)

Does anyone have suggestions on the best way, using Python, to compare a polygon against a list of polygons to identify areas that overlap?


texts: Creating an offset line with the rgeos package

texts:
I discovered the ST_OffsetCurve function in PostGIS, it is based on GEOS.

I'd like to compute the same in the R software, using the package rgeos. I couldn't find a same-named function though.
Has anybody an idea if there is a similar function in rgeos or any other GIS package?


texts: gdal_calc works but I get a Python error at the end of each process?

texts:
I have a problem while using gdal_calc by a .bat file.
Here is the code :

  call "C:\OSGeo4w64\bin\o4w_env.bat"
  gdal_calc -A"D:\tuilage_script_1\tuilage\ndvi\4_810000_6370000_820000_6360000_ndvi.tif" --calc="A" --outfile="D:\tuilage_script_1\tuilage\ndvi\test_gdal_cal3.tif"

It first loads OSGEO environment, then it copies the input raster "A" in an output raster "outfile". The operation is done, but I get the following error messsage : "Python has stopped working", see picture enclosed

I already read this topic :
GDAL_CALC works but I get a python error at the end of each process that prevents automation

But the solution doesn't fit with the aim of my batch file that must be automated to be run from an office basic script for more than 100 operations.

I am using 1.11.2 Gdal version with 2.8 Qgis Wien, 2.7 Python version on Windows-7.


texts: Mapproxy source from WMTS bounding box issue

texts:
I am trying to create an offline version from here:
https://wmts.maptiler.com/aHR0cHM6Ly9ubHMudGlsZXNlcnZlci5jb20vZnBzVVpieHRndGtuL2luZGV4Lmpzb24/wmts
and here is my current config:
services:
  demo:
  wms:
    md:
      title: OS 10k 1888
  tms:
    origin: 'nw'
layers:
  - name: osgb10k1888
    title: OS 10k 1888
    sources: [osgb10k1888_cache]
caches:
  osgb10k1888_cache:
    sources: [osgb10k1888_wmts]
    grids: [webmercator]
sources:
  osgb10k1888_wmts:
    type: tile
    grid: webmercator
    url: https://nls-0.tileserver.com/fpsUZbxtgtkn/%(z)s/%(x)s/%(y)s.jpg
grids:
  webmercator:
    base: GLOBAL_WEBMERCATOR

This all seems to work fine with the resulting levels and resolutions matching the source correct:
Level and Resolutions
Level   Resolution
0   78271.51696402048
1   39135.75848201024
2   19567.87924100512
3   9783.93962050256
4   4891.96981025128
5   2445.98490512564
6   1222.99245256282
7   611.49622628141
8   305.748113140705
9   152.8740565703525
10  76.43702828517625
11  38.21851414258813
12  19.109257071294063
13  9.554628535647032
14  4.777314267823516
15  2.388657133911758
16  1.194328566955879
17  0.5971642834779395
18  0.29858214173896974

however the resulting bounding box is:
-20037508.342789244, -20037508.342789244, 20037508.342789244, 20037508.342789244

I have tried adding the bbox to the grid based on the source:
<ows:LowerCorner>-7.793077 49.852539</ows:LowerCorner>
<ows:UpperCorner>1.790425 60.894042</ows:UpperCorner>

but this the screws up all the resolutions and causes the map to be misaligned
do i need to add all the res into the grid and specify the coords instead of using the web global default?


texts: Showing the current location in ionic with OpenLayers

texts:
I try to show a road in a map + I want the user to be able to see himself on the map.
I have the following
constructor(){
 this.initUserPosition()
}

init(target = 'myMapid', zoom: number) {
private _map: Map
vectorLayer: VectorLayer
view: View
zoom: number

this.view = new View({
  center: [0, 0],
  zoom: zoom,
})

this._map = new Map({
  controls: [],
  target: target,
  layers: [
    new TileLayer({
      source: new OSM(),
    }),
  ],
  view: this.view,
})

new VectorLayer({
  map: this._map,
  source: new VectorSource({
    features: [this.positionFeature],
  }),
})
}
}

I did add the user location like the following.
the onUserPositionChanged$ is an observable that returns my coordinates like this {lng: string, lat: string}
initUserPosition() {
  this.positionFeature = new Feature()
  this.positionFeature.setStyle(
    new Style({
      image: new CircleStyle({
        radius: 6,
        fill: new Fill({
          color: '#3399CC',
        }),
        stroke: new Stroke({
          color: '#fff',
          width: 2,
        }),
      }),
    }),
  )

  this._geolocationService.onUserPositionChanged$
    .pipe(takeUntil(this._unsubscribeAll))
    .subscribe((coordinates) => {
      this.positionFeature.setGeometry(
        coordinates ? new Point([coordinates.lng, coordinates.lat]) : null,
      )
    })
}

Can somebody point me out what I did wrong here?
I already did follow this documentation https://openlayers.org/en/latest/examples/geolocation.html


texts: Calculate distance from GPX in OpenLayers

texts:
I have a little OL script to display GPX Routes in a OSM map. That works so far, even with a little gimmick.

What I miss is a continuous distance. So from trkpt 0 to trkpt 1 + (trkpt 1 to trkpt 2) + (trkpt 2 to trkpt 3) etc.

For the gimmick I use the array keys. Each array key is an (internal) point where data is calculated. In OpenLayers is an array with all coordinates + the elevation. I would like to extend this array by the distance. Whether the array values should be the continuous distance or that from last-trkpt to next-trkpt, I would have to think twice. Maybe a value with last-trkpt to next-trkpt is probably better.

Is there a solution without OL npm. I use only the ol.css and ol.js v5.3.0.

HTML:

<input type="file" id="files" name="files[]" >

Javascript:

var gpxFormat = new ol.format.GPX();
var gpxFeatures;
var coordsarray = [];
var GPXgeometry = [];
function handleFileSelect(evt) {

gpxLayer.getSource().clear();

var files = evt.target.files;

var output = [];
for (var i = 0, f; f = files[i]; i++) {
var reader = new FileReader();
reader.readAsText(files[i], "UTF-8");
reader.onload = function (evt) {
gpxFeatures = gpxFormat.readFeatures(evt.target.result,{
dataProjection:'EPSG:4326',
featureProjection:'EPSG:3857'
});
gpxLayer.getSource().addFeatures(gpxFeatures);

var features = gpxLayer.getSource().getFeatures();
features.forEach(function(feature) {
GPXgeometry = feature.getGeometry();
coordsarray = feature.getGeometry().getCoordinates();
});

//--- ??? Start

var i = 0;
var distance = 0;
coordsarray[0][i].push(distance):

// for ???

var lat0 = coordsarray[0][i][0];
var lon0 = coordsarray[0][i][1];
var latlon0 = [lat0, lon0];
i++;
var lat1 = coordsarray[0][i][0];
var lon1 = coordsarray[0][i][1];
var latlon1 = [lat1, lon1];

// calculate distance
// ???

coordsarray[0][i].push(distance):
i--;

//--- ??? End

}

var gpxLayer = new ol.layer.Vector({
source: new ol.source.Vector({

})
});

var map = new ol.Map({
target: 'map',
layers: [
new ol.layer.Tile({
source: new ol.source.OSM()
}),
gpxLayer
]
});

document.getElementById('files').addEventListener('change', handleFileSelect, false);

Example of var coordsarray:

console.log(coordsarray[0][0][0]); // -37937.014545402875
console.log(coordsarray[0][0][1]); // 5354259.990656373
console.log(coordsarray[0][0][2]); // 188.5
console.log(coordsarray[0][1][0]); // -37936.012669985736
console.log(coordsarray[0][1][1]); // 5354231.092259876
console.log(coordsarray[0][1][2]); // 188.5

EDIT
GPX Example (near equador)

<?xml version="1.0" encoding="UTF-8"?>

<gpx
 xmlns="http://www.topografix.com/GPX/1/1"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
 xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd"
 creator="BRouter-1.4.11" version="1.1">
 <trk>
  <name>brouter_trekking_0</name>
  <trkseg>
   <trkpt lon="11.487967" lat="0.332443"><ele>319.0</ele></trkpt>
   <trkpt lon="11.488969" lat="0.332369"><ele>318.0</ele></trkpt>
   <trkpt lon="11.494446" lat="0.331905"><ele>309.0</ele></trkpt>
   <trkpt lon="11.499806" lat="0.331473"><ele>317.0</ele></trkpt>
  </trkseg>
 </trk>
</gpx>

Javascript:

var startpoint = 0;
var endpoint = startpoint + 2;
var distance = new ol.geom.LineString(coordsarray[0].slice(startpoint,endpoint)).transform('EPSG:4326','EPSG:3857').getLength();
console.log(distance); // 13065083.453663955 ???


texts: How to create a Hillshade model with multicolor "alà Imhof"?

texts:
I would like to make a Hillshade-Model with multicolor "alà Imhof". That means, I dont want a "usual" gray hillshade, I would like to use grey-blue and grey-violett  from different sun directions.

Any idea how I could do that with qgis? Are there any plugins?


texts: Can I extract circular buffers from polygon layers?

texts:
I have been working in QGIS with a large and detailed shapefile where remaining forest fragments are represented as polygons. I would like to extract circlular buffers, of various radii, from this vector layer using R. The forest fragment file and a point I need to buffer can be found here, and a radius I am interested in is 1.128km.

There seem to be ways to extract circular buffers from rasters, but not from polygon layers.  Is the only way for me to perform this operation to first convert my polygon layer into raster form in GRASS, saga-gis, or GDAL?


texts: Join attributes by location - false information

texts:
I've been working on an analysis involving animal vehicle collisions. The data has been passed to r-stats and I keep noticing peculiar things in the cross referencing of the data. I've created highway segments of interest and created a buffer around these. I used the join attributes by location so that I could cross reference the Segment ID onto the Wildlife Collision Data. The Wildlife Collision Data are points that identify the date of a collision and what animal was involved. The highway segments  were created using KDE+ software and they are hotspot segments and coldspot segments. The details are not so important other than to note that each hotspot or coldspot segment has a Segment ID and that is the information I need to be tied to the Wildlife Collision Data.

I've tried several options - different geometric predicates, selecting different fields to add, different joint type (one-to-many v. one-to-one), and so on. I still get errors in the data. The figure below shows the buffered polygon and the points of wildlife collision data contained within.

I know there is a mistake in QGIS because I select the animal column, for example, and note that there are NA's in the in the table of joined output layer. I sort by the segment ID and then select all the data for that segment. If I zoom to that segment I can select the wildlife records to cross compare. The data is inconsistent.

The figure below shows the resulting data table that lists "HORSE" as one of the animals occurring in the segment, but this is FALSE. There's 18 actual wildlife collision records in the segment selected in the figure shown (name = KDE 43 10), but only 14 results total.  There is no HORSE in that segment, but there are records of ELK yet the output doesn't catch these records. Furthermore, some of the joined rows give no output in the view pane. If I select some rows and zoom to selection, I get nothing. There's nothing in these.

I don't know if this is a bug that should be reported and I've never reported a bug. I thought I would post this to see if I'm missing something here.

One thing I've learned is that projection errors are often a likely cause. I decided to look at this data in an older version - QGIS 2.18 - and discovered that the wildlife record layer was not projected correctly. Okay...I fixed the projection. In QGIS 3 it said it had the correct projection and the data fit exactly where it was supposed to fit. I thought that this might be the source of the problem, but the same issue occurs in a re-run of the join attributes. I summarized more detail in the output and you can see in the figure below that the data table in the output records DEER in the appended animal_2 field where it is listed as a HORSE or MOOSE in the animal field. Nothing is making sense.


texts:
I copied a python script from another answered question on this site.  I am using arcgis 10.2.2.  It returns an error and produces an empty shapefile; it has headers but no data.  Can someone look at the script below and identify the problem?

Runtime error
Traceback (most recent call last):
  File "<string>", line 14, in <module>
AttributeError: 'file' object has no attribute 'readLine'

import arcpy

outPath = r"S:\GIS data\Matthew Damage areas" #output path for the feature class
featureClass = "roadSegments.shp" #the name of your shapefile
spatialRef = arcpy.SpatialReference("WGS 1984") #the spatial reference you want; set it to whatever you want

arcpy.CreateFeatureclass_management(outPath, featureClass, "POLYLINE", " ", "DISABLED", "DISABLED", spatialRef) #sets up the shapefile we're using

arcpy.AddField_management(outPath + "\\" + featureClass, "segmentID", "SHORT", 15) #I assume that each segment has some sort of ID number that you might want to import into

roadArray = arcpy.Array() #this will hold the coordinates for each line segment until we write it to the shapefile.

roadFile = open(r"S:\GIS data\Matthew Damage areas\roadlines.csv") #the path to your road lines file
headerLine = roadFile.readLine()
indexList = headerLine.split(",")

#now we need the index locations for the columns that contain the latitude and longitude points for the starting and ending points for each line segment.  I'm also assuming that each segment has some sort of ID value in your segment.  Change the value inside the parenthesis to reflect whatever you named your column for those values.

idValueIndex = indexList.index('segmentID')
fromLatValueIndex = indexList.index('fromLatitude')
fromLongValueIndex = indexList.index('fromLongitude')
toLatValueIndex = indexList.index('toLatitude')
toLongValueIndex = indexList.index('toLongitude')

roadCursor = arcpy.InsertCursor(outPath + "\\" + featureClass)

for fileLine in roadFile.readlines():
    readLine = fileLine.split(",")
    idNumber = int(readLine[idValueIndex])
    fromLat = float(readLine[fromLatValueIndex])
    fromLong = float(readLine[fromLongValueIndex])
    toLat = float(readLine[toLatValueIndex])
    toLong = float(readLine[toLongValueIndex])

    #now let's add the points and draw the lines
    fromVertex = arcpy.Point(fromLong, fromLat)
    roadArray.add(fromVertex)
    toVertex = arcpy.Point(toLong, toLat)
    roadArray.add(toVertex)

    roadCursor = arcpy.InsertCursor(outPath + "\\" + featureClass)
    feature = roadCursor.newRow()
    polyline = arcpy.Polyline(roadArray, spatialRef)
    feature.shape = polyline
    feature.setValue("segmentID", idNumber) #adds the segmentID in the attribute table

    roadCursor.insertRow(feature)
    roadArray.removeAll()

del roadCursor

print "Finished!"
print arcpy.GetMessages()


texts:
In QGIS 3.4.8 I can label vertices with their vertex number by using the following expression for a Font marker symbol:

@geometry_point_num

But, the expression only seems to count the vertices within the current view - so as I zoom in and out the numbers change.

Is there a way to make this work for each whole polyline, irrespective of the view, so the numbers do not change?  So I can, for example, zoom into a line and know I am looking at the fifth vertex along the line.


texts: QGIS server crashes with GetPrint request

texts:
I'm running QGIS server on Ubuntu 14.04 and until now it worked without problems.

Now I want to create print composer output, GetProjectSettings request shows composer template "plantilla".

I request a print output like this:

http://217.126.85.135/cgi-bin/seguimiento_gis/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetPrint&TEMPLATE=plantilla&CRS=EPSG:4326&map0:EXTENT=34.5987,-10.9848,44.5071,4.59801&map0:ROTATION=0&FORMAT=pdf&DPI=300

But fcgid crashes with a coredump:

mod_fcgid: error reading data, FastCGI server closed connection

Update: Using a more simple template works and GetPrint request gives a valid response. Seems to be a problem of texts and HTML in the template "plantilla". Any idea how to debug this error and get to the roots of the problem with the template?


texts: Merge ImageCollections with different pixel size

texts:
I want to combine 2 different Image Collections (i.e. Landsat8 and Sentinel-2), and I want to combine them to a unique ImageCollection (and then run the function qualityMosaic()).
However, they have different pixel size and I think that is giving me the error. I tried resampling but it does not work:

var S2 = ee.ImageCollection('COPERNICUS/S2_SR')

var L8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')

var resample = function (image){
  return image.resample('bilinear').reproject({
    crs: L8.first().projection().crs(),
    scale: 20
})}

var L8 = L8.map(resample)

var S2 = S2.map(resample)

var collection = L8.merge(S2)
Map.addLayer(collection)

Any thoughts?


texts: Creating Demographic Map based on Block Groups with Custom Data using ArcGIS for Desktop?

texts:
I am trying to create a custom demographic map based on block groups that fit the following criteria within a 1 mile radius;

Population = 8,000
Median HH Income = 50k or less
Hispanic Pop = 25% or greater

Is it even possible to select a criteria like this without randomly placing 1 mile radius rings all over the mxd to locate this info?

So far I have created a thematic map using stacked block group layers where the base layer is a graduated color scheme of the income that is less than 50k, then the other two layers for population are dot density. Not quite the look I am going for, but starting at the basic levels for me.


texts: How to import qgis.core in a Osgeo4W installation?

texts:
when I try to import qgis core module from python:
from qgis.core import *
I get this error:
 ImportError: DLL load failed: The operating system cannot run %1.

Searching the web I found that I need to point my PATH to the qgis installation, so I did that:
set PATH=C:\OSGeo4Windows\apps\qgis-dev\;%PATH%

But I still get the error.

obs, i'm using: Osgeo4w Qgis 1.7, Windows 7, Aptana Studio 3 with Pydev.


texts: ArcPy - Select by Attribute for string doesn't work

texts:
I am developing a Python script, in which I want to select features based on an attribute, which is a string. I tried for hours and used a hundred different ways, but I always get the Error 000358 that the query statement is not correct.
I work with a feature class, not shapefile. Assume that my field is called "Fruit" and the value I want to select is "Apple".
What I already tried:
whereClause = """"{} = {}""".format("Fruit", "'Apple'")
whereClause = """"{} = '{}'""".format("Fruit", "Apple")
whereClause = """"{0} = '{1}'""".format("Fruit", "Apple")

fruit = "Fruit"
value = "Apple"
whereClause = """"{} = '{}'""".format(fruit, value)

arcpy.management.SelectLayerByAttribute(myFeatureClass, "NEW_SELECTION", whereClause)

I also did some other Select By Attribute queries, but always using Integer or other numeric fields and I never had problems. Seems that I cannot get it to work with strings.
The error message I get is:
Traceback (most recent call last):
  File "<string>", line 718, in execute
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 8754, in SelectLayerByAttribute
    raise e
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 8751, in SelectLayerByAttribute
    retval = convertArcObjectToPythonObject(gp.SelectLayerByAttribute_management(*gp_fixargs((in_layer_or_view, selection_type, where_clause, invert_where_clause), True)))
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing\_base.py", line 511, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
arcgisscripting.ExecuteError:  ERROR 000358: Invalid expression. Error with executing of SelectLayerByAttribute.


texts: Determining if Point Feature is in More Than One Polygon?

texts: I want to know if a point feature resides inside more than one polygon. I'm using data driven pages and I cannot have the same point feature in multiple frames. There can be over 2,000 point features I'm working with so I need to find a way to locate exactly where a point resides on multiple polygons.


texts: Finding location entering coordinates in QGIS?

texts:
I am new to QGIS.

Can I find a location on the map I'm working with by entering its coordinates?

If yes, where do I enter these coordinates?


texts:
Is it possible to draw one shapeline or one point and repeat it in a specific direction with a specific distance for a specific amount? Like the pattern function in AutoCAD? I'm working with QGIS.
An application example, to digitalize the gouges on a field.


texts: I have uploaded a feature to the PostGIS db and visualize it by mapguide. Now I would like to know how can I update (add) features from database directly. Can I add a new row from the web form directly? Can the mapguide help to generate the new feature? Thanks.


texts: QGIS 3.10 (and newer). Random points inside polygon (Variable) tool is missing?

texts:
Using the current long term release version of QGIS (v3.10), I see the 'Random Points Inside Polygons (variable)' tool is changed. It used to look like this.
The key field here is the 'number field'. Which allows you to select a column in the attribute table of that layer in order to generate a corresponding number of points based on the values of that field.

Now the tool is replaced by simply 'Random Points inside polygons' which looks like this:

The 'Number field' input box is missing.  What happened?  Has the old tool been renamed?


texts: ERROR 000732: Output Location: Dataset results does not exist or is not supported?

texts:
I'm trying to use

  MosaicToNewRaster

but a lot of errors are appears.
One of them is:

  ERROR 000732: Output Location: Dataset results does not exist or is not supported

My code is:

  import arcpy

    from arcpy import env
    env.workspace = "C:\Users\user\Documents\ArcGIS\teste"

    arcpy.MosaicToNewRaster_management("SG-21-Z-D.tif;SG-22-Y-C.tif","C:\Users\user\Documents\ArcGIS","marco.tif","","1_BIT","0.000833333353511989","1")

Should I change something?


texts:
Is there a way to remove shadows of piechart of GeoServer (Eastwood library)? GeoServer always makes shadow in piechart, while the original Google piechart API does not have shadow.

Note: The piecharts in the map above are created with 2 rules, each with a piechart <sld:ExternalGraphic>.


texts: Authentication Error 499 using ArcGIS API for Python?

texts:
Using the ArcGIS API for Python, when I create my GIS object:

gis = GIS(config["agol"]["agol_gis_url"], config["agol"]["agol_gis_user"], config["agol"]["agol_gis_pass"])

I can verify that it connects by checking if I have a token:

token = gis._con.token
print(token)

But If I try to create a FeatureLayer and search, I get an authentication error (499):

p = FeatureLayer(config["agol"]["resource_locs_fs_url"])
result = gis.content.search("id:{0}".format(p.properties.serviceItemId))

Error running result = gis.content.search("id:{0}".format(p.properties.serviceItemId)) in 'C:\test\main.py'. Message: (Error Code: 499)

If I go into ArcGIS Online and share the item to Everyone, the code works fine.  However, the item needs to be secured, so that is not an option.  Shouldn't my gis object grant me access to the item?


texts: Projectionless shapefile with unknown source

texts:
I have a bunch of maps that are output from a land use model in .asc format, but none have a projection file associated. The person who originally ran the LU model is long gone, so I don't have access to the files that were the original input to the model. I've gone through http://spatialreference.org and ad-hoc tried to fit my data to all the projections that looked remotely worth trying, as suggested here and in links within.

Details: Using ArcInfo 10.0. Data is for China. Gauss-Kruger projections came close to looking correct (after on-the-fly projecting on top of my known reference data), but were still several hundred km off in places and somewhat distorted. I think my only recourse is to manually project these projectionless shapefiles - which really sounds like a non-solution to me. Anyone have thoughts on a better fix?

Here's the extent info
Top:    5924000
Left:  -2638000
Right:  2210000
Bottom: 1870000

Thanks


texts: QGIS export degrees with decimal minutes?

texts:
I've been asked to export some coordinates (from WGS84 zone 31N to WGS84) in degrees with decimal minutes.

I've been manipulating an output csv in LibreOffice calc for some while & it seems terribly convoluted (extracting parts of strings & recombining them etc). Perfect breeding ground for errors... Just wondered if anyone had a better way?


texts: How to open a ".sos" file without crashing QGIS?

texts:
I am trying to open a ".sos" file in QGIS. I have QGIS installed using OSGeo4W and added the gdal-sosi package. This should do the trick.

But QGIS just shuts down when I try to open the file. It does not give me any error message.

Do anyone know how I can solve this? I know I can do this in ArcGIS, but I prefer to do it in open source.


texts:
I have 10 fields in my shapefile, and I want to use Python (3.7) to access the value of one field in QGIS.

When I run my code (below) I get a KeyError: '6' (index position of the field), implying that the field lookup worked, but that once I try to access the value, the field doesn't exist.

f = layer.getFeature(fid)
index = layer.fields().lookupField('length')
print("length: ", f[index])

I also tried:

f = layer.getFeature(fid)
index = layer.fields().indexFromName('length')
print("length: ", f[index])

Both throw the same error.

Sample code uses just in a print statement, but I need to use the index number at various points/for different methods.

Clarification: the field look up returns index 6, which is the correct index corresponding to the field length. It's only when I go to use that index number that I then get a Key Error stating that the field doesn't exist.

point_layer = iface.addVectorLayer("...file path.shp", "points", 'ogr')
main_layer = iface.activeLayer()

...more code making lists...

for v_id in list:
    idx = point_layer.fields().lookupField('segment_id')
    v = point_layer.getFeature(v_id)

    segment_id = v[idx]
    s = main_layer.getFeature(segment_id)

    idx = layer.fields().lookupField('length')

    print("length: ", s.attributes()[idx]) #throws IndexError: list index out of range
    print("length: ", s[idx]) #throws KeyError

Image: Works through console - I know the field index is right and the value exists


texts: Zooming in QGIS stops at particular scale

texts:
I seem to have a problem with the zooming levels in QGIS 2.6 Brighton. When I go to zoom in to a vector or raster layer, e.g. for digitizing purposes, the scale seems to be fixed at 1:2,256, I can't zoom in any further than this scale. My previous version of QGIS, 2.2.0 did not have this problem.

Is there a fix for this?


texts: Showing only legend items that are actually displayed on map in ArcGIS Desktop?

texts: Using ArcGIS Desktop, is there a way for the legend, to only show the subitems that are displayed on the map?


texts: how to process modis NDVI data range to normal -1 to 1

texts: I downloaded MODIS Vegetation Indices 16-day L3 Global 500m (MOD13A1) from which i only loaded NDVI value. but the range is 5.25478e+06 to 7.76892e+07. how should I convert it into normal NDVI range of -1 to 1


texts:
I'm using Geoserver 2.6.2. I need to use a SVG file as an external graphic to style a points layer. When I open the WMS layer in QGIS, the points are not displayed. This is not happening if I use a PNG file.

Here is an extract of my SLD code:

<PointSymbolizer>
  <Graphic>
    <ExternalGraphic>
      <OnlineResource xlink:type="simple" xlink:href="http://myserver/img/c_3_3_8.svg"/>
      <Format>image/svg+xml</Format>
    </ExternalGraphic>
    <Size>40</Size>
  </Graphic>
</PointSymbolizer>

I'm wondering why the image format has to be svg+xml.

Is an additional XML file needed?


texts: Breaking up 50mb KML file to under 10mb files?

texts:
I have a large 50mb KML file and we would like the file to be broken into smaller KML files so they can be placed on ArcMap.

How do I do that?


texts: Add Geography Column to an existing table

texts:
I Have an existing table GPSPosition which has a coma seperated LatLong column. I have attempted to add a geography column using the query:

 alter table GPSPositions add geographyColumn as geography::STGeomFromText('POINT('+convert(varchar(20),RIGHT(LatLong, Charindex(',', LatLong) - 1))+' '+convert(varchar(20),LEFT(LatLong, Charindex(',', LatLong) - 1) )+')',4326)

This creates the column but I get problems with the table when I attempt to query it:

  System.FormatException: 24141: A number is expected at position 14 of the >input. The input has ,-1.8635.

How would I construct the query to add the column correctly?


texts: OpenLayers 2 Prevent redraw of specific layer

texts:
In a map with multiple layers I have one WMS-Layer I do not want to redraw on normal map-interaction (zoom, pan). I want to initiate the redraw after a XMLHttpRequest. How can I prevent the redraw for pan/zoom?

I was trying to remove the refresh-listener from the layer, but this has no effect. (I am not sure if this is useful in this context)

Looking into the source I have found the OpenLayers.Map.moveTo() function. Is there a way to configure the layer, to be excepted in the moveTo()?


texts: How to Integrate rivers with forest in ArcGIS for Desktop?

texts:
I have 2 type of inputs, Rivers and Forest. Rivers must not move at all (just freeze), only forest should be forced to attach to rivers within a tolerance of 10-20m.

I want to integrate those layer so they perfectly feat each-other, but when i use the integrate tool in arcgis within a xy tolerance, my rivers moves too in some cases.


texts: Installing QGIS on Mac OS X 10.13.2 with Python 3.4

texts: My first time installing QGIS.  ReadMe says QGIS uses Python 2.7. So, before I start Do I need to remove Python 3.4 or can I install 2.7 as well?


texts: Replicating "Join Attributes by Location" using pyQgis

texts:
I would like to replicate the behaviour of the "Join Attributes by Location" function of the Qgis GUI (Vector >> Data Management Tools >> Join Attributes by Location) using a pyQgis script. I attempt this with the following call:

processing.runalg("qgis:joinattributesbylocation",
                  "path/to/points.shp",
                  "path/to/polygons.shp",
                  u'intersects',
                  0,
                  0,
                  '',
                  1,
                  "path/to/output.shp")

This works except that the resulting output.shp does not contain all the fields of points.shp, it contains only the first two, which happen to be longitude and latitude.

When I use the GUI operation, the resulting output.shp retains all the fields of points.shp and polygons.shp. How do I ensure that all fields of the target vector appear in the output vector?