texts: Cluster GeoJSON (linestrings and markers) into circles using Leaflet

texts:
I have different hiking routes (in geoJSON format) displayed on a map. They contain linestrings, multilinestrings and markers. Currently, it looks like this: http://www.hikethisway.today/Sweden.html.

I would like to cluster the routes in circles where the size of the circle is determined by the length of the route (have a distance attribute in the geoJSON's). Something like this: https://atendesigngroup.com/blog/openlayers-plus
So when zooming out to a certain level, the circle clusters appear and when clicking on a circle cluster the route (line + markers) appear.

I'm currently using Leaflet and the desired output in the second link is using Openlayers to create the visualization. Is it possible to achieve the same result using Leaflet or will I have to move to Openlayers?


texts: Dealing with different MSSQL geometry types in ArcGIS for Desktop?

texts:
I have a Microsoft SQL Server 2008 spatial database holding both point and polygon geometry data types. Displaying the data with the software which captures the data is not a problem, but ArcGIS 10.3.1 cannot display both points and polyons in a single layer.

Is there any way in which I can teach ArcGIS to display both datatypes within a single layer?


texts: How to load more than the standard tools to the grass toolbox?

texts:
I need to know how I get all the grass tools into the grass toolbox of QGIS. Is there an easy way to do this? I tried the sextante toolbox. After an update sextante doesnt work (python errors). For instance: I need the r.series tool.

mfg


texts: Fastest point-in-polygon solution with billions of points over thousands of polygons

texts:
Problem definition is quite simple: I need to achieve point-in-polygon detection for:

billions of lat/lon points
over thousands of polygons

I need to be able, for every point, to determine if it is in one or several of the polygons available (in the thousands). Currently we are using PostGIS with ST_Contains but the whole things needs 3 good entire days (roughly 80 hours) to compute from start to finish.

Is there any solution (free or not) that would offer significant speed improvement for this workload over PostGIS? It can be anything from an analytical database (Exasol is one example that comes to mind, however although hugely faster than Postgres as a database, I'm not sure if it's any faster for geographical computations) to a dedicated C/C++/Rust piece of code.


texts:
I'm working with GDAL using C++.

I've got a code to create a black-and-white image:

GDALAllRegister();
GDALDriver *pDriverTiff = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset *pNewDS;
pNewDS = pDriverTiff->Create(output, nCols, nRows, 1, GDT_UInt16, NULL);
double *transform = geotiff->get_transform();
pNewDS->SetProjection(geotiff->get_proj().c_str());
pNewDS->SetGeoTransform(geotiff->get_transform());
float *newRow = (float *) CPLMalloc(sizeof(float)*nCols);
for (int i = 0; i < nRows; i++){
    for (int j = 0; j < nCols; j++){
        newRow[j] = matrix[i * nCols + j];
    }
    pNewDS->GetRasterBand(1)->RasterIO(GF_Write, 0, i, nCols, 1, newRow, nCols, 1, GDT_Float32, 0, 0);
}
GDALClose(pNewDS);
GDALDestroyDriverManager();

And it works flawlessly.

But now I want to create a colored RGB image. I'm trying to do it this way:

GDALAllRegister();
GDALDriver *pDriverTiff = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset *pNewDS;
char* options[] = {"PHOTOMETRIC=RGB", "PROFILE=GeoTIFF"};
pNewDS = pDriverTiff->Create(output, nCols, nRows, 3, GDT_UInt16, options);
double *transform = geotiff->get_transform();
pNewDS->SetProjection(geotiff->get_proj().c_str());
pNewDS->SetGeoTransform(geotiff->get_transform());
float *newRow = (float *) CPLMalloc(sizeof(float)*nCols);
for (int i = 0; i < nRows; i++){
    for (int j = 0; j < nCols; j++){
        newRow[j] = matrix[i * nCols + j];
    }
    pNewDS->GetRasterBand(1)->RasterIO(GF_Write, 0, i, nCols, 1, newRow, nCols, 1, GDT_Float32, 0, 0);
}
for (int i = 0; i < nRows; i++){
    for (int j = 0; j < nCols; j++){
        newRow[j] = matrix[i * nCols + j];
    }
    pNewDS->GetRasterBand(2)->RasterIO(GF_Write, 0, i, nCols, 1, newRow, nCols, 1, GDT_Float32, 0, 0);
}
for (int i = 0; i < nRows; i++){
    for (int j = 0; j < nCols; j++){
        newRow[j] = matrix[i * nCols + j];
    }
    pNewDS->GetRasterBand(3)->RasterIO(GF_Write, 0, i, nCols, 1, newRow, nCols, 1, GDT_Float32, 0, 0);
}
pNewDS->GetRasterBand(1)->SetColorInterpretation(GCI_RedBand);
pNewDS->GetRasterBand(2)->SetColorInterpretation(GCI_GreenBand);
pNewDS->GetRasterBand(3)->SetColorInterpretation(GCI_BlueBand);
GDALClose(pNewDS);
GDALDestroyDriverManager();

But the image is still black-and-white. How to make it colored?


texts:
I have some countour datasets where some parts were cut by labels, like this:

Can these parts be regenerated? What I've tried is to create a straigth line between the dangling ends, but there are many cases where the created line overlaps other contours.

Another option I've tried is generate contours with the QGIS 'Countour plugin' but to make it work, I must use a high value of point thinning with a result far from what I've expected.


texts: how to fix projection in QGIS from earth data

texts:
After downloading and opening natural earth data in QGIS, the result shows like this

I have tried uncheck simplify geometry but it doesn’t work. I think the problem is related with globe projections but i have never change settings or add any phyton code to do that.

sorry, may be basci but im just started with GIS softhware!


texts: Making a Leaflet line feature more 'grabbable'

texts:
Is there a way to increase the buffer that makes a line clickable in Leaflet?

I know you can increase the width of the line, but the symbology of my feature should not be altered in this case; it should remain exactly as shown in the screenshots below.

For instance, here the line is clickable:

But here it is not, even though the cursor is only a few pixels off:

In OpenLayers I can achieve this pretty easily with the available styling options, which you can tweak to add an invisible buffer, but in Leaflet this does not seem to be possible.

Are there any simple ways to achieve this effect?

I am usinga GeoJSON, in case that makes a difference.


texts: Selecting lat | long | name with given "polygon(,,,,,,,,)" in PostGIS using Python


texts: Raster stacking with naturally ordered layers gives ValueError: need at least one array to concatenate?

texts:
I need to create a multiple bands .tiff file, also known as rasters stack.
I ordered the files to be stacked (from the folder) using natsort and it works.

How to use the method 'os.walk' with the ordered list (ord_list), because
I get this error:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/tmp/tmpAP_F30.py", line 20, in <module>
    stackRast = np.dstack(rasters)
  File "/usr/lib/python2.7/dist-packages/numpy/lib/shape_base.py", line 368, in dstack
    return _nx.concatenate([atleast_3d(_m) for _m in tup], 2)
ValueError: need at least one array to concatenate

I wrote my script doing copy/paste of some stuff in the net, I'm not skilled in Python yet.

import os
import natsort
import sys
import numpy as np
import gdal
import gdalnumeric as gd
from scipy.signal import savgol_filter

rasters = list()

# natural sorting
sourceDir = "/home/giacomo/Desktop/eee/"
unord_list = set(os.listdir(sourceDir))
ord_list = ('\n'.join(natsort.natsorted(unord_list)))

# layer stacking
for folder, subs, files in os.walk(ord_list):
    for filename in files:
        aSrc = gd.LoadFile(os.path.join(folder,filename))
        rasters.append(aSrc)

stackRast = np.dstack(rasters)

How can I avoid getting this error?


texts: GDAL_Translate - ERROR 6: Too many command options - Batch Processing

texts:
I am trying to batch convert a set of JPEGs into ASCII text files.  Below is the code I am running (I added the "..." for privacy, it does not actually appear in the code).

gdal_translate -b 1 -of AAIGrid "C:/Users/.../Cropped Pixel Profile Images/*.jpg" "C:/Users/.../Cropped Pixel Profile Images/ASCII Converted Files/*.txt"

The message I am getting back is ERROR 6: Too many command options 'C:/Users/.../Cropped Pixel Profile Images/-10EV.jpg'.

I would just do this whole process through the "translate batch mode" GUI in QGIS, except that I need to add -b 1 to the code, which can't be done for batch processed items (maybe that's a potential bug/ feature request for future versions of Q).

I am also looking for a way to execute this within QGIS, so that I can run the same command on both OSX and Windows. (If it is platform specific, what is the same command for terminal in OSX).

Anybody have any idea why I am running into this error?


texts: Identifying topological relationships using PostGIS

texts:
I'm trying to perform a basic typology analysis using PostGIS. My objective is to find any polygons that touch other polygons. To do this, I thought the ST_GetFaceEdges would work
(reference). I want to check every polygon in my database and list all the other polygons that touch it. In the image below I would expect the result to say that two of the polygons (which are buildings) touch one building, and the results of the other 4 to say that they touch 0 polygons.

However, I'm having some difficulty understanding what to do. When I tried to copy the example there were a few parts of it that I did not understand.

-- Returns the sequence, edge id
-- , and geometry of the edges that bound face 1
-- If you just need geom and seq, can use ST_GetFaceGeometry
SELECT t.seq, t.edge, geom
FROM topology.ST_GetFaceEdges('tt',1) As t(seq,edge)
    INNER JOIN tt.edge AS e ON abs(t.edge) = e.edge_id;

I'm not sure whether topology is the name of a table, column or a part of the function. I assumed that it was the table, but I'm not sure.

-- try out on sample of data
-- Table is called 'TestArea', column used is 'fid_1'
SELECT t.seq, t.edge, geom
FROM  TestArea.ST_GetFaceEdges('fid_1', 1) As t(seq,edge)
    INNER JOIN tt.edge AS e ON abs(t.edge) = e.edge_id;

I'm also not sure what the function of the inner join is - does this join the result to the original object?


texts: proj4js gives wrong result for EPSG:3375

texts:
Faced with very strange issue. Trying to do coordinate transformation from EPSG:4326 to EPSG:3375 using proj4js and getting wrong result.

This is my code:

var epsg4326 = '+proj=longlat +datum=WGS84 +no_defs'

var epsg3375 = '+proj=omerc +lat_0=4 +lonc=102.25 +alpha=323.0257964666666 +k=0.99984 +x_0=804671 +y_0=0 +no_uoff +gamma=323.1301023611111 +ellps=GRS80 +units=m +no_defs'

var transformation = proj4(epsg4326, epsg3375)
var input = [101.70979078430528, 3.06268465621428];
var output = transformation.forward(input);
// 744626.6428051005, -103618.62269607143
// expected: 412597.532715, 338944.957259

This is a JSFiddle workspace for live demo of the issue.

I have created an issue at the proj4js github repo.

Any ideas on that?

Could you please suggest alternative JS libraries I can use to do this transformation in browser? It looks like many other libraries depends on proj4js which is not perfect as I found out.


texts: Create isochrones from streets with distance values?

texts:
I have made isochrone lines using GRASS and QGIS, as shown. The lines show various distances to fire stations on a street network.

Now I would like to turn these lines into polygons that look like isolines using QGIS. The lines to polygon tool does not work as expected. Is there a way?


texts:
I am using ArcGlobe.
How can I add an image that is just on the side, and not associated with the globe?
I want to add a color scale, but I don't know how.


texts: Why can't I connect to my PostGIS database?

texts:
It's the second time that I install PostgreSQL with the PostGIS extension and it runs as long as I don't restart the PC. After I have restarted the PC QGIS and PgAdmin4 can't connect to the database. How can I fix that?

I have PostgreSQL 9.6.1-1 and PostGIS 2.3.1 and Windows 7.


texts: Visualizing shapefile in CAD the same as GIS

texts:
I have a shapefile of conifer trees in QGIS 3.16 that I exported as dxf and opened in Global Mapper. In QGIS, the area covered by conifers is easily recongnizable with the blue color as in IMG1.

In dxf format, only the lines around the areas are visible, so it is not possibile to distinguish easily whether the tree polygons are inside or outside of these lines as in IMG2.

The issue is that when I select an item in the dxf file, it selects all the area of that polygon and not only the separate blue areas as in IMG3. As shown, the white areas inside the big blue polygon are also selected in dxf, while they are no conifers there.

Is there a way to visualize the polygons in dxf format as they are in GIS? Can I change something in my polygon so it appears the same as I see it in GIS for my boss who sees it in CAD?


texts: I had 2 versions of pandas (22 and 23) in my Anaconda environment for Python 3.6. One of the pandas looks like it was installed through Python 2.7 pip. Basically, I am unable to import pandas at all into my code. It has been recommended that I completely uninstall and reinstall Anaconda. This seems like overkill, but I am going to try it. Nothing I attempted using conda or the Anaconda environment console seemed to work. My question is, how do I safely uninstall and reinstall Anaconda without messing up my ArcGIS Pro environment or my Python 3.6 (do I need to uninstall and reinstall that too?!) Help!


texts: Calculating midpoint between every two closest points using ArcGIS Pro

texts:
I am looking to try and identify the midpoint between two closest points in my layer in ArcGIS Pro.
The reason being that I work in the utilities industry and I am trying to calculate the angle between every structure in the network. Using the Near tool I was able to calculate it, but there are instances where structures are side by side. For this instance I would look at the two very close structures as one structure for the purpose of finding the angle in between. The structures side by side give me an angle of 0, 90 or 180 when calculating the angle via the near tool.
I have subsetted them to only focus on those "problem structures" and tried things such as aggregating points, buffering and intersection etc but am now stuck.


texts:
I was wondering if someone could point me in the direction of display the output of multiple gps devices at the same time in QGIS.

I have two GPS devices that are both outputting NMEA sentences through either a com port or over a TCP connection.

I can get either one to plot easily but not both together.


texts: rioxarray and vector data (including geocube)

texts:
I currently got some problem with rioxarray and geocube. I hope, this is the right place to ask.
I am using a vector file for the German state of Brandenburg, so it has a massive hole in it. The file only has this single shape in it. When I create a raster file with:
    bb_raster = make_geocube(vector_data=bb_frame, resolution=(-60, 60), fill=-999)

But when I plot the data the plot is empty (white). But np.all(bb_raster['id'].data == -999) returns False
Because of its size it is much larger, than a tile from the sentinel satellites for example.
When I want to mask it with a sentinel-2 raster with:
scl = scl.rio.clip(bb_frame)

the

ValueError: No valid geometry objects found for rasterize

is returned instead. But still it is "allowed" to ...rio.reproject() a sentinel-2 raster with the rasterized shape.
But when I plot these reprojected files, then all the plots are empty (white). Although valid data seems to be included the data xarray.DataArray.
An example:
    xarray.DataArray

        band: 1y: 4078x: 3894

        array([[[241, 241, 241, ...,   4,   4,   4],
                [241, 241, 241, ...,   4,   4,   4],
                [241, 241, 241, ...,   4,   4,   4],
                ...,
                [241, 241, 241, ..., 241, 241, 241],
                [241, 241, 241, ..., 241, 241, 241],
                [241, 241, 241, ..., 241, 241, 241]]], dtype=uint8)

        Coordinates:
            x
            (x)
            float64
            2.502e+05 2.502e+05 ... 4.838e+05
            y
            (y)
            float64
            5.935e+06 5.935e+06 ... 5.69e+06
            band
            (band)
            int32
            1
            spatial_ref
            ()
            int32
            0
        Attributes:

        scale_factor :
            1.0
        add_offset :
            0.0
        grid_mapping :
            spatial_ref
        _FillValue :
            241

What might be the cause and the solution to this weired behavior?


texts: Connecting multiple cities on a map using MapQuest route API

texts:
I am working on any optimal vehicle routing problem.
I want to start map routing on the road starting from
one city A
to city B to City C
 to city D
back to city A

Based on lat/lng of these cities

I have the cities already worked out just want to display it on map-quest map?

This is sound as it is fundamental what mapping solution will offer out of the box but I am had googled this for days now with no particular solution

Can somebody point out the right direction to go?

I am using MapQuest.


texts: Using ESRI.ArcGIS.Display.ITemplate to draw cutom lines?

texts:
I am using ESRI.ArcGIS.Display.ITemplate to draw custom lines such as dot/dash lines on a map control.

Then, I used below codes

ESRI.ArcGIS.Display.ITemplate templateCls = new ESRI.ArcGIS.Display.TemplateClass();
templateCls.Interval = 1;
templateCls.AddPatternElement(1, 3);

The below image 1 will show the result for this.

But, I am want to the line include round elements (not rectangle) as the below image 2.

Are there anyone know how to do that ?

Thank you very much !


texts:
I'm a complete beginner trying to update some fields in an attribute table (i.e. an ArcMap Shapefile) using Python tool. My goal it to update the values from this column to another date notation.

Right now the fields look like this:

  91/02/25
   93/01/21

  96/06/14
  ...

The dates need to be replaced to the following format MM/DD/YYY. There are also some spaces (as shown above in the example) that I would like to replace to something like "No date available".

This is what I wrote:

for row in updateCursor:

    theWrongDate = row.RPT_DTE
    theRightDateList = theWrongDate.split("/")
    year = '19' + str(theRightDateList[0])
    month = theRightDateList[1]
    day = theRightDateList[2]
    theWrongDate = month + '/' + day + '/' + year
    updateCursor.updateRow(row)

and Python gives me an error message:

"month = theRightDateList[1]
IndexError: list index out of range"


texts: Problem reprojecting data from personal geodatabase using QGIS

texts:
I have a dataset in the form of a personal geodatabase (mdb) and am interested in using certain feature classes for further processing.

My problem is - I am using USGS data orthoimagery which comes in UTM coordinate system and am unable to reproject my geodatabase feature class (layer) to UTM. The same task is easily done in ArcGIS, but I want to use QGIS as there are other features of QGIS that I intend to use later on.

What I have tried:
1) Save geodatabase layer (feature class) as shp and try to reproject that shp --> the reproject task does not proceed beyond 2% completion and I end up with a 1KB shp file

2) Reproject mdb layer (Feature class) directly to UTM --> It does not allow me to do this either.

3) Create shapefile from mdb layer (feature class). Define projection for the shp file. Reproject this shp file. --> this also results similar to (1)

So, for some reason, projecting layers that have their origin from an ESRI mdb do not seem to be reprojected in QGIS.

Any help/tips on this would be welcome.

Thanks in advance.


texts: Are there asset management systems for PostGIS?

texts: Are there asset management systems that use PostGIS for asset storage?


texts:
I’m trying to detrend a DEM in QGIS.
The best guidance that I could find on the internet was this question : Removing elevation trend over sloped surfaces?

I couldn’t understand that so I would like to know if someone could explain how to remove the tilt of a DEM.
I know that the best approach is to fit a plane and then subtract this from the original DEM. But I can’t understand the criteria to create this plane and how I perform the subtraction.
I can use python if needed.


texts: Using gdal_grid tool - set a defined extent

texts: I want to use the gdal_grid interpolation tool. Is it possible to set a defined extent?


texts: Using Yelp API lat/long in Folium -- CRS needed?

texts: I understand the default crs for Folium is EPSG3857. I am pulling lat/long data from a yelp query and am not sure if I need to be converting to a GeoSeries and assigning a CRS to the GeoSeries. It seems to be drawing fine, but I want to make sure I'm doing everything in the correct way. Does anyone have any ideas?


texts: Splitting and converting a large tif files into multiple .hgt files

texts:
I'm quite new to GIS format and conversions, so be patient with a newbie. I got stuck on a file conversion and I'm looking for some help this to fix my scripts. I have a big 5 GB tif file with a precision of 5 meters and these features:

Dimensions X: 27176 Y:.49586 Band: 1
Origin 555121,4.4447e+06
Pixel dimension 5,-5
No value:-9999
AREA_OR_POINT=Area
STATISTICS_MAXIMUM=2237.2060546875
STATISTICS_MEAN=520.43425935058
STATISTICS_MINIMUM=-2.8499999046326
STATISTICS_STDDEV=429.22835999081
Data type: Float32
+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs
Layer extension (original layer projection)
555120.6040000000502914,4196769.7539999997243285 : 691000.6040000000502914,4444699.7539999997243285

I'm getting that each hgt file covers 1 degree lat and 1 degree lon, and has to have a size of 1201x1201.

My intent is to convert this single tif into one or more hgt files. I suppose that given the tif dimensions (X: 27176 Y:.49586) I have first to split into multiple tif files and then convert into multiple hgt file. To do this I'm using GDAL scripts from QGIS:

If I try a direct conversion to hgt I get this warning and error:

Warning 1: The corner coordinates of the source are not properly aligned on plain latitude/longitude boundaries.
ERROR 1: Image dimensions should be 1201x1201, 3601x3601 or 1801x3601.

So I defined these scripts below:

Step 1:

gdalwarp -of GTiff -srcnodata 32767 -t_srs "+proj=utm +zone=33 +datum=WGS84" -rcs -order 3 -ts 1201 1201 -multi C:/Users/C012426/Desktop/geocoding/5m_33N/DEM_5m_33N.tif C:/Users/C012426/Desktop/geocoding/5m_33N/DEM_5m_33N_warped.tif

Step 2:

gdal_translate -of SRTMHGT C:/Users/C012426/Desktop/geocoding/5m_33N/DEM_5m_33N_warped.tif C:/Users/C012426/Desktop/geocoding/5m_33N/DEM_5m_33N_warped.hgt

What I get instead is one file only named: N4196770E55.hgt

while I was expecting to get 27176 / 2101 = about 22 hgt files...

What is wrong in my theory and most of all in my scripts ?


texts: Raster projection is inverted Python-GDAL

texts:
I am carrying out a program that, based on the data from weather stations, performs an interpolation and generates a georeferenced raster with these data.

The program works well interpolating and saving the raster, however, the projection of the raster is inverted.

This is the image of how the raster is stored, the green dots are the location of the weather stations, which should coincide with the vertices of the raster image.

I know that interpolation works well because when I use Matplotlib to visualize the image it looks like this:

Because of this, the problem must be when I do the raster projection that corresponds to this part of the code:

#Saving raster
xmin,ymin,xmax,ymax = [xi.min(),yi.min(),xi.max(),yi.max()]
nrows,ncols = np.shape(zi)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)

output_raster = gdal.GetDriverByName('GTiff').Create('myraster.tif',ncols, nrows, 1 ,gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)  # Specify its coordinates
srs = osr.SpatialReference()                 # Establish its coordinate encoding
srs.ImportFromEPSG(4326)                     # This one specifies WGS84 lat long.

output_raster.SetProjection( srs.ExportToWkt() )   # Exports the coordinate system
                                               # to the file
output_raster.GetRasterBand(1).WriteArray(zi)   # Writes my array to the raster

output_raster.FlushCache()

output_raster= None

I've been modifying the program for a long time and I can't find the problem, I hope you can help me.


texts: OpenGeo Suite 4.8 - service unavailable after enable CORS edits and restart (Jetty)

texts:
I have Open Geo Suite 4.8 installed on a server.  I tried to enable CORS by following the guide given here. I have the jetty-servlets file in the directory given in the Boundless instructions.

After restarting  OpenGeo Jetty (version 7.6.13), I now get a 503- Service Unavailable message when I try to view the map services in our web app.

Is there another step that I am missing? The steps given in the Boundless page are pretty straighforward.

I already reviewed the answer given to Enabling CORS in GeoServer (jetty)?, but I am still getting a service unavailable error.


texts: GRASS: import PostGIS vector layer into SQLite database

texts:
I am working with vector layers (points) composed of more than 100 000 geometries. These layers are originally in a PostGis database and I import them with v.in.ogr, using a formulation like:

v.in.ogr \
    input="PG:host=127.0.0.1 dbname=gisdb user=user password=secret" \
    layer=schema.points \
    output=points \
    -o --overwrite

This command is successful, however, operations with the resulting points layer are remarkably slow. For instance, this layer must be later loaded in R and that operation alone is taking over 3 hours.

After some back-and-forth testing, I identified the problem: v.in.ogr is not saving its output in the internal GRASS SQLite database. Knowing this, I can create a copy of the original output, that is automatically inserted into the SQLite database, e.g.:

g.copy vector=points,points_mirror

Importing the points_mirror layer to R takes 3 minutes instead of 3 hours, for instance.

While mirroring vectors imported with v.in.ogr solves the issue, it is an inelegant solution that also takes up unnecessary space. I would prefer to instruct v.in.ogr somehow to force it to use the SQLite database. I could not find any reference in the manual, is it possible to do?


texts: Spatial Query between polyline features and point features within specified Buffer

texts:
I'm new to GIS (so I'm lost on this) and looking for a C# code sample to help me do 2 types of spatial queries. I have 2 Layers only. First layer is pipelines, which are polyline shapes so I can select multiple pipeline segments to make up an entire pipeline. Second Layer is facilities, which are point shapes.

Query 1: Select Pipeline Segments in the Pipeline Layer (Polylines) and find the facilities (in the Facility Layer) that are within a certain distance of selected Pipeline Segments.

Query 2: Select a single Facility in the Facility Layer (Point) and find the Pipeline Segments (Polylines) that are within a certain distance of the Facility Point.

I can do these queries via the ArcMap interface. I select the Source Features in the Attribute Table, then Selection Menu, Select by Location, select the "other" Layer as the Target Layer and for Spatial Selection Method I select "Target layer(s) features are within a distance of the Source layer feature". But now I need to do this in C# code for my ArcMap AddIn.


texts:
I am using OpenStreetMap for showing my user a map within my web frontend. A user can mark on the map his parcel.

I would like to automatically get from the marked area information about the height difference on the parcel. Is this possible with OpenStreetMap?


texts: Error: basic_string::_S_create, (trying to use .ecw files)

texts: I have a problem, after resolving my problem with: .ecw is not a supported data file. I am in a new one. When I try to load an .ecw  file in QGIS it says: basic_string::_S_create


texts: Limiting available values in dropdown field of QGIS form according to value selected in previous field?

texts:
I'm building a rather complex entry form for a layer stored in a PostGIS database. I am following this usage scenario:

On several occasions I have dropdown fields which lets the user select a type of artefact (populated through value relation to another table in a database), followed immediately with the field letting them select a subtype (also populated through value relation to another table in a database). I would like a choice of subtypes in the dropdown to be limited by the previous selection of a type.

For example, if type "ceramics" is selected the subtype dropdown needs to be limited to only options with a key value from 1 to 10.

If type "metals" is selected the subtype dropdown needs to be limited to only options with a key value from 11 to 20.

And so on.

I feel relatively fluent in SQL and qgis expression builder, and i would guess this could be achieved trough filter expression option on value relation widget, but I wouldn't know how to invoke the values from a form which is currently being filled.

I suspect some python code is needed? - I'm not fluent in python but I guess I could adapt an example of code to my needs by renaming the fields etc.

Question is similar to Field Value Mapping QGIS but the answers there didnt point me out in directions I could apply with my current knowledge.


texts: Exporting Network Database (from Network Analyst Tool) to where it would still be built on another computer from download

texts: I have a Network Database built within GIS for a transit provider with the stops, lines, etc. all built onto a routable network. I would like to export the routable network within the Network Database to another computer but when I exported it, it only had the layers (not the Network Database that is routable / works within Network Analyst).


texts:
I am trying to perform some scale-tests with some data on GEE. I have the following code where I use the MODIS-NDVI dataset at 500m. What I want to do, is to perform the same analysis but on 25000m x 25000m. Any ideas how to change the pixel size of MODIS-NDVI dataset?

var modisNDVI = ee.ImageCollection ('MODIS/MCD43A4_006_NDVI');
var modiscollection = ee.ImageCollection (modisNDVI.filterDate('2002-03-01','2004-10-31'));

    var chart = ui.Chart.image.doySeriesByYear(
        modiscollection, 'NDVI', region, ee.Reducer.mean(), 500);

    print(chart);

    var clip = modiscollection.mean().clip(region);
    Map.addLayer (clip, {min:0.0,max:1,palette:['FFFFFF','CC9966','CC9900','996600','33CC00','009900','006600','000000']},'NDVI');


texts:
I've run into this problem with Leaflet Map. I am trying to loop over some addresses, make a call for each loop to get the lat-long, plot on mymap, and group the markers to bound them using fitBounds
createRequest(projectUrl, getAddresses, handleErrors, projectInit)
    .then( (jsonAddresses) =>{
        // jsonAddresses = [{building_num: 100, street_name: 'main', ... }, {}, ...]

        let state_province = ''
        let markers        = []
        let projectsGroup  = L.featureGroup().addTo(mymap)   // Creating an empty group to add markers to later on

        for (let address of jsonAddresses )
        {

            let theAddress  = address.building_num + ' ' +
                              address.street_direction + ' ' +
                              address.street_name + ' ' +
                              address.street_type + ' ' +
                              address.civic_address_city + ' '+
                              address.state_province + ' ' +
                              address.country

            addressInit = {'address': theAddress ,
                           'buidlingNum': address.building_num,
                           'streetName': address.street_name,
                           'city': address.civic_address_city
                          }

            // Using geoApify API to get lat-long for each address
            geoCodingUrl = `${geoapifyURL}${theAddress}&apiKey=${GEOAPIFY_API_KEY}`

            // createRequest makes asynchronous call and returns a json response
            createRequest(geoCodingUrl, getLatLongGeoApify, handleErrors, addressInit )
            .then( (response) => {
                latLongAddress = response
                latLongAddress.address = address  // address comes from the for-loop

                latLongs.push(latLongAddress)   // latLongAddress is a list of objects [{<lat>,<lon>}, {},...] returned by createRequest
                return latLongs
            })
            .then( (result)=> {
                let marker = L.marker( [latitude,longitude] ).addTo(mymap)
                marker.addTo(projectsGroup)

                markers.push(marker)

                for (let project of result)  // project represents the lat-long of the project, e.g. project = [<lat>,<lon>]
                {
                    marker.setLatLng( [project.lat, project.lon] ).bindPopup("Try!") // Setting the marker coordinates
                    marker.addTo(projectsGroup)   // Adding the marker to the projectGroup
                }
                mymap.setZoom(12)

            })

        }

        // test_marker is to create another marker that cannot fit in the
        // default view, which necessitates fitBounds to work properly
        const test_marker = L.marker( [48.354495, -122.0905781] ).addTo(mymap)
        test_marker.addTo(projectsGroup)

        mymap.fitBounds(projectsGroup.getBounds() )

        return latLongs

    })

});

When I look at my markers in the console, it looks like an empty list, but when I open it, there are 4 markers in it.  I know this must have to do with asynchronous nature of the call and the fact that the console is live, so initially it does not have any data and then it receives the data.  But, I cannot figure out how I can add markers to the featureGroup in the async-proper way so I can bound them


texts:
I am trying to create a median or average stack profile on a DEM in QGIS but I can't find any tool or plugin for it. This feature is implemented in GMT 1).
Can anybody give me a hint of how it can be done in QGIS?


texts:
I am running caching on a server using ArcGIS Server 10. It seems like it is holding SDE connections, older than 2 weeks, to the point that it has maximized connections. I am not sure why it is not freeing up resources.

What I can I do to make sure that it releases SDE connections.

Thanks
Jay


texts: Overlay shapefile on multiband raster using ggRGB and ggplot2 in R

texts:
I am trying to use ggplot2 to map a polygon feature class over a raster brick (created from multiple Landsat band images). I'm able to map the polygons and the raster separately, but when I try to have them together (following the last code chunk at http://bleutner.github.io/RStoolbox/rstbx-docu/ggRGB.html RStudio crashes (even after restarting the computer).

ggRGB attempt:

fallowFields <- rgdal::readOGR(dsn = fgdb, layer = 'Fallowing2018')
fallowFieldsWGS84 <- spTransform(fallowFields, crs(L8_20180420_brick))

fallowFieldsDataFrame <- fallowFieldsWGS %>% fortify

landsatAprilCropped <- raster::crop(L8_20180420_brick, raster::extent(fallowFields))

p <- ggplot()
p + ggRGB(img = landsatAprilCropped,
          r = 3,
          g = 2,
          b = 1,
          stretch = 'hist',
          ggLayer = TRUE) +
geom_polygon(col = 'Black',
             fill = 'grey70',
             data = fallowFieldsDataFrame,
             aes(x = long, y = lat, group = group))

replacing the data frame (fallowFieldsDataFrame) with the shapefile (fallowFieldsWGS84) works, but it warps the Landsat image.

The cropped image and shapefile are in the same projection, but the shapefile seems to be warping the raster. Is there a way to plot the shapefile over the raster without warping it?

It doesn't get warped with plotRGB, but I was hoping to use ggRGB so that I could make a more elegant map with gglot2.

plotRGB(landsatAprilCropped,
        r = 3, g = 2, b = 1,
        stretch = 'hist', axes = TRUE,
        main = 'Fallow Fields Overlaid on Landsat Image: 4/20/2018')

plot(bardFallowFields, col = 'transparent', border = 'red', add = TRUE)


texts: OpenLayers 3: projection for Yandex.Maps

texts:
I need 3 layers: OSM, Google and Yandex. My view, OSM and Google layer in EPSG:3857 projection, but Yandex layer seems to be in EPSG:3395, so it was displayed not correctly. I've defined the new projection via proj4js and the layer has moved significantly, but it still doesn't match the other layers.

My code is:

proj4.defs('EPSG:3395', '+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs');

map = new ol.Map({
  target: 'olmap',
  layers: [
    new ol.layer.Tile({
      title: 'Yandex',
      source: new ol.source.XYZ({
        url: '...',
        projection: 'EPSG:3395'
      })
    }),
    new ol.layer.Tile({
      title: 'Google',
      source: new ol.source.XYZ({
        url: '...'
      })
    }),
    new ol.layer.Tile({
      title: 'OSM',
      source: new ol.source.XYZ({
        url: '...'
      })
    })
  ],
  view: new ol.View({
    projection: 'EPSG:3857',
    extent: ol.proj.transformExtent([-180, -80, 180, 80], 'EPSG:4326', 'EPSG:3857')
  })
});

The definition of EPSG:3395 I took from here: https://epsg.io/3395 .

OSM layer:

Yandex.Maps layer:

Where is my mistake?


texts: Why do labels in exports using HERE tiles appear tiny?

texts:
I am trying to export images from the QGIS print composer, the exported image has a basemap that is loaded from an XML file following this questions method. The export doesn't look at all like the image I see in the print composer see below and I'm trying to understand why.

What I'm trying to do is simply make a decent looking map without obvious tiling with readable labels. Perhaps I have to use a satellite only image and remake with my own labels. The map is a HERE nokia hybrid map

This is a jpeg at 150 dpi, note the tiny labels in comparison to the screenshot at the bottom
.

This is a png export at 300 dpi, which is standard for the print composer, note the labels, and tiles.

Finally, this is a screenshot of the map as it looks in the print composer with the labels readable at this zoom scale.


texts: What is difference between simultaneous users and simultaneous connections in ArcGIS for Server 10.1?

texts:
What is the difference between simultaneous users and simultaneous connections.  We're looking at Workgroup standard vs Enterprise Standard.  Here's what we want to do:  Use Arc Gis Server to serve data to app. 30 vehicles using Ipads via Sierra Wireless 4g modems.  This data would be maps for Right of Way, outage locations etc.  We would also like to be able for users to send data to other users.

Thanks,

Russell


texts: Can animation exports be longer than 10 seconds?

texts:
I'm exporting an animation to .avi and they all turn out to be 10 seconds long. Is there a way to change this? To have the animation play for a longer duration?

ArcGIS 10.2.2


texts: Locus Maps and Rmaps Android maps support an offline sqlite Format called Rmaps. I am trying to find the format details and how its organized. I couldn't, is there any dependable documentation online?


texts: Deleted features not synchronizing in one-way (ArcSDE to File Geodatabase) geodatabase replication?

texts:
I am testing three one way  (ArcSDE 10 -> File Geodatabase) replicas.  Each replica is derived from a layer created by selection by location from a citywide road dataset, so that each of the three layers is distinct from one another.

When I add features (all line layers) to the ArcSDE parent replica, all the features appear in the children after synchronization, as I would expect.  The problem occurs when I delete any new test feature that is beyond the extent of the original dataset.  If I first add a new feature beyond the dataset extent to the parent, this feature will be added to the child, but then when I delete the same feature in the parent, this feature will remain in the child after synchronization.

I do not have this issue when testing replicas that use the entire (non-selected) dataset.  I have tried some different combinations of extent settings in the replica creation and in the environment settings but none have worked yet.  I suspect this is related to the fact that these test features do not intersect the original extent...but I can add features outside of the extent, and I do not have this problem when testing  a replica using the 'full' dataset.

Any ideas?


texts:
I cached my map hoping to use it with openlayers. This is my code below.

    var layers = [
    new ol.layer.Image({
      source: new ol.source.ImageWMS({
        url: 'http://0.0.0.0:8080/geoserver/wms',
        params: {'LAYERS': 'silapha:basemap'},
        ratio: 1,
        serverType: 'geoserver'
      })
    })
  ];

  var map = new ol.Map({
    layers: layers,
    target: 'map',
    view: new ol.View({
      center: [0, 869099],
      zoom: 7
    })
  });

This still renders the old uncached layer. I enable direct wms integration  as well but I can't just render the cached map.

This this is the url to the cached map and the preview works well.
http://0.0.0.0:8080/geoserver/gwc/demo/silapha:basemap?gridSet=EPSG:900913&format=image/png

What have i done wrong, what should add to the code to get the cache tile rendering instead of the uncached one.


texts: Add polygon attributes to points fails to create layer

texts:
I'm trying to add polygon attributes to a point file with a lot (~370,000) of points in QGIS v3.0.0.  I first tried using the point sampling tool, which I've used with success pulling the values from raster files, but this became unresponsive and remained so for several hours.

A search pointed me to the Saga add polygon attributes to points from the processing toolbox, and reported it to be much faster than join attributes by location see this post by @Underdark link.  This completed quickly but failed to create or add a layer.  I got a similar error as @Tali in this question.  I followed the multipart to singlepart suggestion, but to no effect- my polygons were all initially single part features.  I also ran v.build.check to see if I had any other topological errors in my polygon file, but found none.

The error message I get is:

I can do other processing tasks with the point layer, and tried both in and out of edit mode.  The tool works fine with the same polygon layer and a much smaller point files (~400 pts), but still failed when I used a ~200k point subset.  I've tried re-running with the output set to save as both as shapefile and as geopackage, but get the same result.
I converted both my point and region files to geopackages and tried again, with the same error. I also tried in a blank map document with only the two files and get the same error.

Any ideas on how to resolve the error with the add polygon attributes to points, or an alternate tool that could accomplish the same thing in a reasonable amount of time?


texts: Port ArcGIS Engine to ArcGIS Runtime SDK for WPF?

texts:
We have a WPF application, which has arcengine controls embedded in it. We are looking at options to replace ArcGIS Engine with someother component which is 64 bit and we are evaluating ArcGIS Runtime SDK for WPF (possibly .net as well).

We do not need an ArcGIS for Server or ArcGIS Online based solution, we dont use web maps or tiles services. We have file geodatabases and need to create, edit and modify feature classes within the gdb.

1**. Can I use ArcGIS Runtime SDKs to build a non-cloud / non - server WPF Solution.i.e. The map control in the WPF sdk will need to access only a local fgdb (not hosted as a service).**


texts:
When I try to run QNEAT3's "OD Matrix as Lines/CSV" function, I am getting a very strange error.
Traceback (most recent call last):
File "C:/Users/dicki/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QNEAT3\algs\OdMatrixFromPointsAsLines.py", line 207, in processAlgorithm
net = Qneat3Network(network, points, strategy, directionFieldName, forwardValue, backwardValue, bothValue, defaultDirection, analysisCrs, speedFieldName, defaultSpeed, tolerance, feedback)
File "C:/Users/dicki/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QNEAT3\Qneat3Framework.py", line 130, in __init__
self.list_tiedPoints = self.director.makeGraph(self.builder, self.list_input_points, self.feedback)
Exception: unknown

Execution failed after 53.70 seconds

Any thoughts? I'm at a loss - I've posted this in a few other places as well.


texts:
I would like to graph my data from a shapefile. The data I wish to display is the sum of the nominal data per each year that it occurred.

I have enclosed a screenshot to illustrate my idea. There is the field of data (string) I wish to graph. There are only A's and W's in the field (125 in total) and I want the W's totaled per year and the A's totaled per year as multiple bars like shown in the graph below. Then I can display those bars on a yearly basis with their respective totals. If that makes sense, please help me. Thank you.


texts: making choropleth maps using PostGIS Windows

texts:
i read an article on chicago tribune blog about making maps,
it consist on "making choropleth maps using PostGIS, TileMill, Mapnik and Google Maps" on Mac environnement, my question is how we can make same thing on windows evironnement.
i want to make a Side-by-side map comparisons like this

Thx


texts:
I would like to perform snapping using Arcpy for a lot of files in all the lines of the shapefiles that don't connect to the nearest other lines, so they do.

I tried to do it in geopandas using Python but I can't make it work.

Since I don't know much about ArcPy I will show you what I have done in geopandas:

import geopandas as gpd
from shapely.ops import snap
s=gpd.read_file(r'C:\Users\user\Desktop\ag_eus')
result=snap(s,s,0.5)

How this code would be converted in ArcPy?

I have done this process in regular Arcmap and worked great but how to do it in ArcPy so it can work with multiple shapefiles in a loop. If you can explain what you do too, would be great.


texts: How to attach timestamp to single track (not points) in GPX file

texts:
I'm trying to build a Torque visualization with CartoDB based on a GPX file I received from a third-party. Unfortunately, it doesn't seem to be formatted quite right as the trkpts don't have timestamps but the trk as a whole does.

Example:

<trk>
    <name>02bb8f71fa569898748f1983d78a2c21</name>
    <desc>fd367fa83bb00443b2614da748fa9a40</desc>
    <timestamp>2013-11-25 18:47</timestamp>
    <trkseg>
      <trkpt lat="28.543573333333335" lon="-81.35272833333333"/>
      <trkpt lat="28.543573333333335" lon="-81.35272833333333"/>
      <trkpt lat="28.54366" lon="-81.35315833333334"/>
      <trkpt lat="28.54366" lon="-81.35315833333334"/>
      <trkpt lat="28.54366" lon="-81.35313666666667"/>
    </trkseg>
</trk>

It should look like this:

<trk>
    <name>02bb8f71fa569898748f1983d78a2c21</name>
    <desc>fd367fa83bb00443b2614da748fa9a40</desc>
    <trkseg>
         <trkpt lat="28.542491666666667" lon="-81.37588166666667">
             <time>2009-10-17T18:37:26Z</time>
         </trkpt>
         <trkpt lat="28.54224333333333" lon="-81.37438">
             <time>2009-10-18T18:37:26Z</time>
         </trkpt>
    </trkseg>
</trk>

Is there a way I can still do a Torque visualization without having times for the individual points? After uploading the data into CartoDB, two tables are produced -- one with track information, one with point information. The track table doesn't even include the timestamp associated with it in the GPX, and Torque isn't an option in the Map wizard.

Here is a screen shot of a row from the tracks table:

And here are two screen shots (table is too wide to fit in one) of the points table.


texts: Error while creating large TIFF with gdal_merge

texts:
I need to merge two TIFFs into one, but the output is larger than 4 GB. So, it needs to be a BigTIFF. But the "BIGTIFF=YES" creation option seems to be ignored. Here's the error:

$ gdal_merge.py -separate -co "BIGTIFF=YES" -co "COMPRESS=DEFLATE" -o merged.tif vv.tif vh.tif
0...10...20...30...40...50ERROR 1: ZIPDecode:ZLib error:
ERROR 1: TIFFReadEncodedStrip() failed.
ERROR 1: merged.tif, band 1: An error occurred while writing a dirty block from GDALRasterBand::RasterIO
...60...70...80...90...100 - done.
ERROR 1: TIFFAppendToStrip:Maximum TIFF file size exceeded. Use BIGTIFF=YES creation option.
ERROR 1: ZIPDecode:ZLib error:
ERROR 1: TIFFReadEncodedStrip() failed.

The output TIFF (which is corrupted) is 4.3 GB whereas the two TIFFs I want to merge are both 1.5 GB.

I also tried without the compression option, but the result is the same. Any idea where this comes from?


texts: How does osmosis take so much disk space?

texts:
I've been trying to use osmosis for importing osm raw data to a db. First, it failed due to lack of sufficient disk space. I made some more space and launched it one more time. That time there was another error and the execution was aborted anyway BUT a huge amount of space was still taken from my hard drive... I had 11 GB before executing an osmosis procedure and now I have like less than 1 GB and it's getting difficult to use the computer anymore. The question is: what files  did osmosis create on my disk during the procedure execution and how can I find them and delete them in order to have my disk space back?

I'm using the following osmosis command for import:

"YourPathToOsmosis\bin\osmosis" --read-pbf file="D:\YourPathToOSMFile\osm_file.osm.pbf" --write-pgsql host="localhost" database="yourDBname" user="postgres" password="yourPW"


texts: Leaflet R - plotting raster with colour palette

texts:
I am trying to plot a raster image using Leaflet, however I am having issues getting the colours right.
I've read in the raster like this and plot the the resulting raster variable.
library(raster)
library(leaflet)
ras<-raster("https://raw.githubusercontent.com/jodyphelan/jodyphelan.github.io/master/img/test.tif")
plot(ras)

But when I try to plot the raster using leaflet and the addRasterImage() function there are two things I notice. Firstly, it looks as its it is doing some kind of antialising between the colours. Secondly, it is putting some colours in that shouldn't be there.
pal<-colorNumeric(ras@legend@colortable,0:255,na.color = "transparent")
leaflet() %>% addRasterImage(ras,colors = pal)

This is the raster produced using plot():

An this is the one produced using leaflet:

The questions are:

Is it possible to disable the antialising-like behaviour?
Why is it changing colours of some squares?

I'm also not entirely sure whether I am creating the pallete in the right way. Any tips on these questions?


texts: Selecting features based on great number of different values in the Attribute Table

texts:
I have many land parcels related to many different cities. I would like to select land parcels related to some of these cities only (about 20 of them).
Is using LIKE/ILIKE + OR like below the only way to do this or is there a faster and more elegant way to do it when dealing?
'ILIKE'-way:
"city" ILIKE 'Paris'
or
"city" ILIKE 'Rome'
or
"city" ILIKE 'Berlin'
or
"city" ILIKE 'Moscow'

Maybe in the form? But just writing 'Paris' OR 'Rome' or 'Berlin' or 'Moscow' in the "city" field doesn't seem to work.


texts: Double rotation in SLD style in Geoserver

texts:
I am facing a problem with Rotation in SLD styling in Geoserver. I have a symbol that consists of two triangles. One of them rotated 180 degrees so that they point towards each other.

So far so good - that works fine for me

But I would also like the total symbol (i.e. the two triangles together) to be rotated based on an attribute value called "Symbolretning"

Whatever I do I cannot get this "dobbeltrotation" to work.

Do any of you a way to solve this? I think the main issue of my problem is where in the SLD code I put the second rotation statement.

Here is my code so far - without the attribute based rotation:


texts:
I am starting to learn to develop stand alone QGIS applications from a book called 'Building Mapping Applications with QGIS'.

Here's the python code stated in the book (lex.py):

import os, os.path, sys

from qgis.core import *
from qgis.gui import *
from PyQt4.QtGui import *
from PyQt4.QtCore import *

class MapExplorer(QMainWindow):
    def __init__(self):
        QMainWindow.__init__(self)
        self.setWindowTitle("Landmark Explorer")
        self.resize(800, 400)

def main():
    QgsApplication.setPrefixPath(os.environ['QGIS_PREFIX'], True)
    QgsApplication.initQgis()

    app = QApplication(sys.argv)

    window = MapExplorer()
    window.show()
    window.raise_()
    app.exec_()
    app.deleteLater()
    QgsApplication.exitQgis()

if __name__ == "__main__":
    main()

Batch file (run.bat):

SET OSGEO4W_ROOT=C:\OSGeo4W64
SET QGIS_PREFIX=%OSGEO4W_ROOT%\apps\qgis
SET PATH=%QGIS_PREFIX%\bin;%OSGEO4W_ROOT%\bin;%PATH%
SET PYTHONPATH=%QGIS_PREFIX%\python;%PYTHONPATH%
SET PYTHONHOME=%OSGEO4W_ROOT%\apps\Python27
python lex.py

However, when i run the batch file. it just displays this:


texts:
In Arcmap I have made a new raster with ArcToolBox -> Spatial Analyst Tools -> Map Algebra -> raster calculator.

Now, months later, I have forgotten the source formula expression.

How can I display?


texts:
What does "Precision" do within "Join Attributes by Location"?

Also, what's the difference between the "Geometric Predicates"?


texts: Draw a spatial grid with the lat/longs at the center of the grid using R

texts:
I looked through multiple posts on the same, and most of them end up drawing a grid connnecting the points or a grid connecting the bounding box of the coordinates.

I need a grid with each of these points at the centre of the grid. This is what I have so far:

library(sf)
library(ggplot2)

points <- data.frame(long = c(76.75,77.25,76.75,77.25,77.75),
                     lat = c(11.25,10.75,10.25,10.25,10.25)) %>% st_as_sf(coords=c('long','lat'), crs=4326)

a <- st_make_grid(points, what="polygons", cellsize=0.5)

ggplot() +
  geom_sf(data = points, color = 'red', size = 1.7) +
  geom_sf(data = a, fill = 'transparent', lwd = 0.3) +
  coord_sf(datum = NA)  +
  labs(x = "") +
  labs(y = "")


texts: Converting point data into a weighted raster

texts:
I am completely re-writing my question as i think it was quite confusing before:

I have a number of different points, each point has a -z value.
I need to create some sort of raster or a series of polygons which buffers these points.. but i have a few complicated rules!

If the adjacent points are of the same z value then the buffer must merge into one polygon.
if the adjacent point is a different z value - then the buffers must be weighted in accordance to the depth of the z values, i.e. if one z-value is -10, and the other is -5, then the buffer of the -10 point will take up 75% of the distance between the 2 points. But because there are many points, and they are all related, this buffering will be different on different sides of the point, see diagram!
The idea of the weighted voronoi diagram has been suggested to me and that would be perfect if it weren't for all of the points being related to all of the points around its!
I know how i would do it on paper - i would divide the space between all points by 10, and then based on the difference between the two points choose where the buffer would start (if there is no difference between the points then the buffer starts in the middle).
Although i cannot work out the formula for determining where the buffer starts.. i know it must be derived from the DIFFERENCE between the two points - but i do not have the mathmatical capacity to work out how i would use this to work out the distance!

I am using Arcmap 9.3.1 and am open to python scripting!
If you need me to try and explain anything further i will try! i am very confused myself!


texts: How to change GK-coordinates to UTM with the field calculator in QGIS?

texts:
i have  a qgis project in ETRS89 (UTM) with some vectorshapes. i already transformes them from GK into UTM. but now the attribute table still show me the GK-coordinates. HOW can i change it? with the field calculator?!

Please help me :)
Lena


texts: QGIS expression builder include relation attributes

texts:
I have a set of symbols based on a table in postgis such as this

id|type|action_1|action_2|geom

In QGIS I have this table symbolised using a rule expression which works effectively.

However this table was setup as a bit of a rush job, as it doesn't adhere to 3NF structure and is inflexible (i.e. what if more than 2 actions needed). As such the structure is now:

Table 1

id|type

Table 2

id|tbl_1_id|action

Where I have a 1-n relationship between table 1 and table 2. I have incorporated the forms and entering data works perfect. However I am trying to incorporate the same rule based symbology as I had previously.

I have tried using the  reffuntion plugin, which doesn't do this (or I haven't been using correctly). I have tried creating a view in my database and joining it to table one so I have a field which I can then symbolise on. However when I edit the data, the view table does update - however the joined table does not update - even after pressing the refresh data button.

Any thoughts on how to symbolise based on two tables OR refreshing data in joined tables would be great.


texts: Get the whole path and name of selected layer from QcomboBox in QGIS plugin

texts:
I am creating a temporary plugin for QGIS using the QGIS plugin builder plugin. Aim is to learn how to make plugins in QGIS.
I have added two comboBoxes that allows the user to select a vector layer and a raster layer respectively. The plugin will process the data specified by the user in these comboboxes. For that I need to extract the entire path and name of the layers selected by the user. This might probably seem quite a basic post but I tried looking up for this, and failed to find a solution even after playing around with several codes. Below is the code for the comboBoxes:

    # Getting the loaded raster layers
    layerlist_rst = ['<Select a raster layer>']
    layers_rst = QgsProject.instance().mapLayers().values()
    for layer in layers_rst:
        if layer.type() == QgsMapLayer.RasterLayer:
            layerlist_rst.append(layer.name())

    # Getting the loaded vector layers
    layerlist_vect = ['<Select a vector layer>']
    layers_vect = QgsProject.instance().mapLayers().values()
    for layer in layers_vect:
        if layer.type() == QgsMapLayer.VectorLayer:
            layerlist_vect.append(layer.name())

    #Clear contents of the comboBox from previous runs
    self.dlg.comboBox_rst.clear() #This combobox loads raster layers
    #Populate the comboBox with names of all the loaded raster layers
    self.dlg.comboBox_rst.addItems(layerlist_rst)
    #Clear contents of the comboBox from previous runs
    self.dlg.comboBox_vect.clear() #This combobox loads the vector layers
    self.dlg.comboBox_vect.addItems(layerlist_vect)

Now I need to get the names and path of the selected layers as text.
Example: If the layer added by user is 'C\Users\XYZ\Rstr.tif'
I want to assign this path to a variable that can be used later in the plugin like:
raster = r'C\Users\XYZ\Rstr.tif'


texts: Save/Export specific contour levels for a KDE (package 'ks' in R) as ESRI shapefile

texts:
I'm performing a kernel (KDE) home range analysis of an individual with the package 'ks' in R. I'd like to avoid using Geospatial Modelling Environment (GME) for implementing the output in ArcGIS, because I'm using QGIS which is not compatible with GME. I'd like run the analysis in R (see example below) and import the output as a shapefile into QGIS.

gps <- read_csv("data/gps.csv")
gps<-data.frame(cbind(gps$x,gps$y))
colnames(gps)<-c('x','y')
gps_hpi<-Hpi(x=gps,pilot="samse",pre="scale")
gps_kde<-kde(x=gps,h=gps_hpi,compute.cont=TRUE)
contourLevels(gps_kde,cont=c(50,75,95))
contourSizes(gps_kde,cont=c(50,75,95),approx=TRUE)
plot(gps_kde, display="filled.contour", cont=c(50,75,95))

These contours are exactly what I want to export as a shapfile, but there is still no reasonable possibility described in the internet to do this.

Adapting @Jeffrey Evans' Formula resulted in a KDE-plot which differed significantly from the original plot.kde function.

r<-raster(extent(440000,465000,5655000,5680000),nrows=nrow(gps_kde$estimate), ncols=ncol(gps_kde$estimate))
r[]<-gps_kde$estimate
r<-flip(r,direction='y')
rcont<-rasterToContour(r,levels=gps_kde$cont)
(cont.values<-gps_kde$cont[grep(paste(c("50","75","95"),collapse="|"),names(gps_kde$cont))])
rcont.gt50<-rcont[which(rcont@data[,1] %in% cont.values),]
par(mfrow=c(1,1))
plot(gps_kde,display="image", main="KDE")
plot(r,main="raster of KDE")
plot(rcont,main="All contours (1-99%)")
plot(rcont.gt50, main="50%, 75% and 95% volume contours")

Further analyzes are not possible with this result. Also the method of @Juan Antonio Roldán Díaz resulted in non-applicable results.

hts<-contourLevels(gps_kde, prob = c(0.5, 0.75, 0.95))
c50<-contourLines(gps_kde$eval.points[[1]],gps_kde$eval.points[[2]],gps_kde$estimate,level=hts[1])
c75<-contourLines(gps_kde$eval.points[[1]],gps_kde$eval.points[[2]],gps_kde$estimate,level=hts[2])
c95<-contourLines(gps_kde$eval.points[[1]],gps_kde$eval.points[[2]],gps_kde$estimate,level=hts[3])
Polyc50<-Polygon(c50[[1]][-1])
Polysc05<-Polygons(list(Polyc50),"c50")
Polyc75<-Polygon(c75[[1]][-1])
Polysc75<-Polygons(list(Polyc75),"c75")
Polyc95<-Polygon(c95[[1]][-1])
Polysc95<-Polygons(list(Polyc95),"c95")
spObj<-SpatialPolygons(list(Polysc05,Polysc75,Polysc95),1:3)
axu.df<-data.frame(ID=1:length(spObj))
rownames(axu.df)<-c("c50","c75","c95")
spDFObj<-SpatialPolygonsDataFrame(spObj, axu.df)
plot(spDFObj)

I have tried many ways (mainly tranforming into raster and then creating contour lines) but they all failed because their results were always significantly different to the original KDE-contours-plot (plot.kde).

How do I export the contour lines of my KDE (50%, 75% and 95%) as a polygon or polyline shapefile so I can load it into QGIS for further analysis (same as shown above)?


texts: How to show data defined centroid fill marker font in the legend?

texts:
I'm using QGIS 3.18.0 on Windows 10 (from OSGeo4W Network Installer).
I have a polygon-shapefile with an attribute table including the attributes "signatur", "r", "g" and "b":

I'm using Simple Fill together with Centriod Fill with Font Marker to distinguish between the polygons:

I'm using Data defined override to set the color of the Font marker to be based on the values in the "r", "g" and "b" fields:

I'm using Data defined override to set the Font Marker character(s) to be based on letters in the "signatur" field:

Looks like this in the QGIS canvas:

Looks the same in the Layout. In the Item Properties - Legend Item Properties  - Symbol settings has the same setup as in Layer Styling. However the Legend patch is either missing the actual Font Marker character (if I have chosen a blank character) or showing an "A" (A seems to be the default selected character) in the default color red:

I want the patches in the legend to show the same Font Marker (characters = "signatur"-value) as in the map, and not like this:

Is that possible?


texts: Finding all distinct values from single column of an attribute table with GDAL

texts:
I am working with GDAL and modify ESRI shapefiles with it. I am searching a possibility how to find all distince values from a single column of an attribute table with GDAL (Python). In SQL it would be SELECT DISTINCT field FROM table. I found how to do that in ArcPy but I prefer to use open source software.

Furthermore, I was thinking about iterating through the column (which are different classes) and calculate their area. Anybody an idea?

for i in range(layer.GetFeatureCount()):
  feature = layer.GetFeature(i)
  #here somehow iterate/sort through "feature.GetField("classes")"


texts: TIN interpolation from LiDAR point data excluding no-data areas (for water)

texts:
I am working with LiDAR data for a specific marsh area. The data is in x,y,z point files. When loading this layer, it is clear that in the creeks of the marsh, no data points were measured with the LiDAR (because of the presence of water in the creeks at time). When I perform a TIN interpolation however, the no data parts (creeks) are still interpolated, which is probably skewing the resulting DTM.
I've tried using breaklines in the TIN interpolation to exclude the creek areas with no data points, but then the resulting image still includes this creek area (with very diverging elevations inside the creek area).
This is what I've done in the TIN interpolation using breaklines
QGIS code revision: 654e76b3cb
Qt version: 5.11.2
GDAL version: 3.1.4
GEOS version: 3.8.1-CAPI-1.13.3
PROJ version: Rel. 6.3.2, May 1st, 2020
Processing algorithm…
Algorithm 'TIN interpolation' starting…
Input parameters:
{ 'EXTENT' : '148996.216800000,149565.618900000,210533.096900000,211123.716200000 [EPSG:31370]', 'INTERPOLATION_DATA' : 'C:/Users/LIDAR/LiDAR_new/BW_2014_clipped_xyz.gpkg|layername=BW_2014_clipped_xyz::~::1::~::-1::~::0::|::MultiLineString?crs=EPSG:31370&field=fid:long(0,0)&uid={c35efefe-93ff-4ee4-a305-18620c0578cf}::~::0::~::0::~::2', 'METHOD' : 0, 'OUTPUT' : 'TEMPORARY_OUTPUT', 'PIXEL_SIZE' : 0.1 }

Execution completed in 23.85 seconds
Results:
{'OUTPUT': 'C:/Users//OUTPUT.tif',
'TRIANGULATION': ''}

Loading resulting layers
Algorithm 'TIN interpolation' finished

Is there a way to do this interpolation without taking into account the no data points, so the resulting DTM is the marsh area, with just an empty/no data area of where the creek is?
(I'm using QGIS 3.16 and have no experience with Python)


texts: Iteration algorithm over geometric attributes of shapefile with arcpyscripting

texts:
I am trying to create a separate skyline shapefile for every feature ID I have in my point file by splitting them into multiple files. The code works fine, except splitting them into multiple shapefiles takes extremely long time. I have 65000 points in the shapefile. I have attached a photo of how I split my files. Once I split them, I point to each shapefile using the below code.

Splitting so many points currently takes about 25-30 hours on my computer. I use QGIS for this because I do not have advanced license for arcpy.Split_analysis tool to work with ArcMap.

Is there a better way than splitting the point file every time?

import arcpy, glob
# reading all the shapefiles
arcpy.env.workspace = "C:\\centrelinepoints"
files = glob.glob("C:\\centrelinepoints\*.shp")

#running the skyline function on each point file separately.

for file in files:

    ## creating a unique name extension for each skyline shape output file based on its feature ID.
    name = file
    parts = name.split('_')
    newname = parts[5]
    final = newname.split('.')
    FID = final[0]

    ## definiing the output layer location and name based on their feature ID.
    output_layer = "H:\\skylineshapes\Skyline_" + FID + ".shp"
    output_angle = "H:\\skylineangles\Skyline_" + FID + ".dbf"

    arcpy.Skyline_3d(file, output_layer, inSurface)

    arcpy.SkylineGraph_3d(file, output_layer, "0", "ADDITIONAL_FIELDS", output_angle)


texts: How do I use QGIS to merge maximum values of multiple TIF files into one "maximum value' raster?

texts:
I have 14 rasters that represent water storage for the same area over a range of years, 2002 to 2015. I am looking for a way to find the maximum value at each pixel location.

I am not looking for the maximum value for each raster; I am looking for the  maximum value for each pixel of 14 rasters at each location.


texts: ArcGIS Pro failed to connect to PostgreSQL 11 database

texts:
I'd like to connect from ArcGIS Pro 2.4.2 to a PostGIS database in PostgreSQL 11.2.
I am able to access via pgAdmin 4 to the respective databse. The link and access with my credentials is working for other systems e.g. QGIS etc. The steps were "insert" "connections" "new database connections" filling in my credentials then I get following error message:

  failed to connect to specified DBMS instance, a database connection
  exception has occurred. underlying DBMS error(FATAL:no pg_hba_conf
  entry fsqlstate=)

As I am able to connect from other systems I don't know what is missing in the pg_hba.conf file.


texts: for loop did not generate a 'BASIN' for each 'CROSSING'

texts:
My code is intended to delineate watersheds from a specified "CROSSING". The purpose is to obtain the total drainage area for a point on a stream/river labeled "CROSSING". My issue is that only one "BASIN" shapefile is generated, containing all "CROSSING" points, at the location furthest downstream.
How can I alter my code to have individual "BASIN" polygons drawn, starting from the "CROSSING" point at the highest elevation?
The image is of the generated shapefile from my current code.
    import arcpy
    from arcpy import env
    from arcpy.sa import *
    #DEFINE PARAMETERS & zLimit
    env.workspace = arcpy.GetParameterAsText(0)     #type is folder
    inSurfaceRasters = arcpy.GetParameterAsText(1)  #type is Raster Layer, multiple
    CROSSINGS = arcpy.GetParameterAsText(2)         #type is Table View
    EXTENT = arcpy.GetParameterAsText(3)            #type is double
    PARCEL = arcpy.GetParameterAsText(4)            #type is Table View
    zLimit = 3.2808
    #CHECK OUT REQUIRED EXTENSION
    arcpy.CheckOutExtension("Spatial")
    #CLIPPING RASTER LAYER TO SMALLER AREA
    inSurfaceRaster = arcpy.MosaicToNewRaster_management(inSurfaceRasters)
    AOI = arcpy.Buffer_analysis(PARCEL,'AOI',10000)
    USR = ExtractByMask(inSurfaceRaster,AOI)
    #GENERATING REQUIRED RASTER LAYERS
    outFill = Fill(USR,zLimit)
    outFlowDir = FlowDirection(outFill)
    outFlowAcc = FlowAccumulation(outFlowDir)
    arcpy.AddXY_management(CROSSINGS)
    #DELINEATION OF WATERSHEDS
    #---------------------------------------------------------- begin loop
    fc = CROSSINGS
    field = "POINT_X"
    cursor = arcpy.SearchCursor(fc)
    for field in cursor:
        outSnapPoint = SnapPourPoint(fc, outFlowAcc, EXTENT)
        watersheds = Watershed(outFlowDir,outSnapPoint,"VALUE")
        arcpy.RasterToPolygon_conversion(watersheds,"BASINS")
        #DELETE EXCESS SHAPEFILES
        arcpy.Delete_management(outSnapPoint)
        arcpy.Delete_management(watersheds)
    #---------------------------------------------------------- end loop
    #DELETE EXCESS SHAPEFILES
    arcpy.Delete_management(AOI)
    arcpy.Delete_management(USR)
    arcpy.Delete_management(outFill)
    arcpy.Delete_management(outFlowDir)
    arcpy.Delete_management(outFlowAcc)
    #LET USER KNOW LAYER IS READY
    print "Layer(s) is(are) ready to be added from "+env.workspace


texts: Tool validation to filter out OID and geometry fields?

texts:
I have a python script tool (arcpy) in a standard toolbox that edits data in fields, but I want some code in tool validation that would make sure that when the user gets to the field parameter, the check box list of fields will never contain OID or geometry.

How might that look in the tool validation code? I have almost no experience tool validation code by the way. I have no code because I don't know what this would look like.


texts:
I am trying to render a map provided by a tile map service based on {x}, {y} and {z} using leaflet.js. So far, I managed to render the map, but it is isn't shown correctly. The tiles are repeated oddly, and the Lat and Lng aren't accurate.

fiddle

For example, Lat: -2.6000285, Lng: 118.015776 with a zoom level of 5 should get the map to show the country Indonesia (see coords at google maps)

But it is not the case in my jsfiddle. Instead, Indonesia can be found around
Lat: 66, Lng: 331, Lat: 66, Lng: -29, and probably at other points every 360 increment/decrement as well. The point is, the lat-lng is not accurate to real world map.

How do I fix this?


texts: How to update the attribute table in layer from dictionaries using python

texts:
I have dictionary with the values:

dict = {0L: {'val1': 60, 'val2': 199, 'val3': 7.719},1L: {'val1': 60, 'val2': 199, 'val3': 7.719}}

I want to update these values in the attribute table of a layer. How can i do that?


texts:
I have OpenLayers 3.0 map with vector layer(Heritage Layer). I am trying to popup the details of the Heritage while clicking on the location of the heritage.
The Problem is: clicking on the symbol of the heritage does not click on the layer.
 var view = new ol.View({
projection: projection,
center: [84.101, 28.44],
zoom: 6.8

});
map = new ol.Map({
layers: layers,
target: 'map',
overlays: [overlay],
controls: ol.control.defaults().extend([
  new ol.control.ScaleLine({
    units: 'metric',
    target: 'scale'
  }),
  zoomToExtentControl,
  new ol.control.FullScreen()
]),
view: view

});
var popup = new ol.Overlay({
element: element,
positioning: 'bottom-center',
stopEvent: false

});
map.addOverlay(popup);

var HeritageLayer = new ol.layer.Tile({
id: 'heritage',
title: 'heritage',
visible: false,
legendPath: 'webatlas:district_wgs_update',
source: new ol.source.TileWMS({
    url: 'http://localhost:8080/geoserver/webatlas/wms?service=WMS',
    hidpi: false,
    params: {
        'VERSION': '1.1.0',
        'LAYERS': 'webatlas:heritage',
    },
    serverType: 'geoserver',
    crossOriginKeyword: 'anonymous'
})
});

map.on('singleclick', function(evt) {
  console.log("hey");

// Attempt to find a feature in one of the visible vector layers
var feature = map.forEachFeatureAtPixel(evt.pixel, function(feature, HeritageLayer) {
    return feature;
});

if (feature) {
  alert("I am an alert box!");
  }

else {
        console.log("Feature Not displayed");

    }

});

Here in this code Feature Not displayed is only displayed while clicking on the heritage layer.


texts:
I am trying to merge all road data from the OS Open Map Local file, to create a single line shapefile of all the roads in the UK. I will then use this data to create huff catchments. I would like to do this all in one step as opposed to spRbind 50 + times. I was able to batch load the shapes into my R environment but am now stuck. Any advice?

shps <- dir("C:/R_Code/MSc_Engagement_Model/1.Data/OS_Road_Local/Roads_Only", "*.shp")

for (shp in shps) assign(shp, readOGR(shp))

roads <- spRbind(HP_Road.shp, HT_Road.shp) #??


texts:
We have recently switched from ArcGIS to QGIS and I'm trying to simply limit the amount of data that shows to a defined area.

In ArcGIS this is simply 'data frame clip'.

Can I do this in QGIS?

I don't want to delete any data, just limit what is shown by clipping to a rectangle or view.


texts: Convert raster_stats list output to array / csv file

texts:
I usually program in R, but recently I have started to use python for most of my spatial programming. My python skills are still limited, and please excuse the simplicity of the question. I am struggling to convert a list, which is an output of a spatial function, into an array / csv file. In particular, I struggle how to disaggregate the list output that i receive from the spatial function.

Here is the spatial function I have used:

from rasterstats import raster_stats
import time

os.chdir("C:/Users/Documents/veg_dt/data")

start = time.time()
stats_n_w = raster_stats("counties.shp", "veg_ndv.tif",stats=['count'], categorical=True)
end = time.time()
print end - start

The function raster_stats counts the occurrence of a categorical variable in a set of polygons. The output of the function is a list, and this is the actual output:

{0: 31728, 1: 35, 2: 31, 3: 30, 4: 29, 5: 40, 6: 44, 7: 82, 8: 85, 9: 142, 10: 20, 11: 58, 12: 76, 'count': 32400, '__fid__': 39989}
{0: 32204, 1: 6, 2: 13, 3: 5, 4: 15, 5: 27, 6: 33, 7: 5, 8: 28, 9: 28, 11: 3, 12: 33, 'count': 32400, '__fid__': 39990}
{0: 31983, 1: 38, 2: 11, 3: 21, 4: 42, 5: 17, 6: 19, 7: 52, 8: 66, 9: 34, 10: 9, 11: 87, 12: 21, 'count': 32400, '__fid__': 39991}
{0: 32400, 'count': 32400, '__fid__': 39992}

Every element of the list counts the categorical variable in a given polygon. The first row shows 31,728 zeros, 35 ones, 31 twos, 30 threes, 29 fours, 40 fives, etc in polygon 39989. The second element has 32,204 zeros, six ones, etc in polygon 39990. I would like to convert the list into a data frame / array / eventually csv file that states the fid (id of the polygon cells) in the first column, then lists the occurrences of zero in the second column, lists the occurrences of ones in the third column, etc. I want to have every element of the list on one row.

Also, to complicate things further, the list is around 60,000 elements and I will do the same operation approximately 50 times for different sets of polygon and tiff images.

Can anyone help with this?


texts: How to set a defined distance between nodes to digitize in qgis?

texts:
Using QGIS (by now I have the 2.4 version) I would like to set some parameter (I don't know if such exists) that would impose a definite or maximal distance between nodes when digitizing a polyline or polygon.

I need this to respect a minimal grain in the mapped features which I have to bring out.

For example I would like to set a digitizing distance between nodes of -let say- 15 m, that is, when digitizing, I will be forced to make a node either all 15 meters (not more, not less) or at least all 15 m (i.e. within few to max. 15 m).

I searched in QGIS menu : Settings > Options > Digitizing, and also in plugins but could not find anything that seems appropriate to me. But maybe I overlooked something ?

Thanks a lot for any help.

NB : I don't want to calculate nor find out distances ! I want to set them.


texts: GeoWebCache and OpenLayers horizontal resolution error

texts:
I have a mapping application (the Arne G. Brekke Bygdebok Collection) which uses OpenLayers to load map tiles from a GeoWebCache server operated by the Statens Kartverk (Norwegian Mapping Authority).

It's worked nicely for two years, but this morning the tiles all load as pink, with errors like "400: Requested horizontal resolution: 0.0439453125 , best match: 0.03333333333333333 exceeds 10% threshold. Perhaps the client is configured with an incorrect set of scales (resolutions), or the DPI setting is off compared to the one in GWC?" if you view the image directly.

Here is a test case I worked out demonstrating the problem:

<!DOCTYPE html>

    <meta charset="UTF-8" />
    <title>Arne G. Brekke Bygdebok Collection</title>

Arne G. Brekke Bygdebok Collection

That's the simplest test case I can work out for it.  Based on that, OpenLayers generates URLs like this one (broken onto multiple lines for clarity):

http://opencache.statkart.no/gatekeeper/gk/gk.open
    ?LAYERS=topo2
    &FORMAT=image%2Fjpeg
    &SERVICE=WMS
    &VERSION=1.1.1
    &REQUEST=GetMap
    &STYLES=
    &SRS=EPSG%3A4326
    &BBOX=11.25,56.25,22.5,67.5
    &WIDTH=256
    &HEIGHT=256

How can I get it to work again?  Bear in mind that I have zero control over the map server -- I can only adjust the JavaScript on my end.


texts:
I have a point layer (stored in SQL Server 2008) that is being served up by a map service in ArcGIS Server 10.1.  I am consuming this service in a web client.

Lets say the points are boats out at sea.
This is mock data, but I want to simulate these boats moving around in a logical pattern.

At the moment, I have an arcpy script that updates the X,Y within a bounding box, and also updates two fields for speed and direction (which the layer is symbolised by)

This is all done using the random module from Python.

It works fine, but does not look realistic.

I have formulated some ideas about having a LUT, with X,Y, Speed, Direction, for each 'boat' that could be referenced to update by.

However, this must have been done before, but I have yet to find any examples.  I would prefer not to spend development time to try to re-invent the wheel, so was wondering if anyone had a similar script for achieving this that they could share as an answer for all?

Reason I choose boats, as its not essential for them to be following a path such as a road.
The simulation does not need to be 99% accurate on how the (insert moving something) actually moves around in the real world, but does simulate movement.

Perhaps if the answer was using a pre-cooked LUT, then the paths could be defined beforehand?


texts: Setting initial extent in Configurable Map Viewer (CMV) via ArcGIS API for JavaScript API?

texts:
I've been playing around with the Configurable Map Viewer (CMV) and would like to set a custom initial extent.

I was thinking I could replace

center:[-96.59179687497497,39.09596293629694]

with

extent: ({xmin:-123.45,ymin:123.45,xmax:123.45,ymax:123.45,spatialReference:{"wkid":102100}})

in viewer.js but either I'm missing something entirely or I'm not doing something quite right.

Could anyone provide an example of setting a custom initial extent in CMV?


texts: I have brought in multiple TIFF (raster) files into QGIS and want to combine them into one file.  I tried this by going to Raster> Misc.> Merge.  When doing this the output file comes out as gray and black and I want to keep the colors of my TIFF files.


texts: Writing fields from point shapefile to text file using ArcPy?

texts:
I have the following script, trying to take a cities shapefile and create an exported text file that will list each feature in the text file with (I can't figure out where I'm going wrong):

CityName1, CoordX, CoordY

CityName2, CoordX, CoordY

import arcpy
arcpy.env.workspace = "C:\\Data"
citiesFile = "C:\\Data\\Cities.shp"
fields = ["SHAPE@XY"]
output = open("C:\\Data\\CO_Cities.csv", "w")
with open("C:\\Data\\CO_Cities.csv", "w") as COCitiesFile:
    with arcpy.da.SearchCursor(citiesFile, "SHAPE@") as cursor:
        point = row[0].getPart(0)
        citiesFile.write("{}, {}\n".format(point.X, point.Y))


texts: How can I convert an Integer field to a Real field in a shapefile?

texts:
How can I convert an Integer field to a Real field in a shapefile ?

I'm searching a tool, or a way to do that (may be under QGis ?!).

Regards


texts: Acquiring administrative country border data with a given resolution

texts:
I want to be able to find out which country a given (lat,long) pair falls into, but I don't need this to be extremely accurate and especially coasts can be quite crude as false positive hits will be fine.

OSM maps with administrative borders such as the set from https://mapzen.com/data/borders/ come very close to what I want, but their resolution is much higher than I need so the files are a magnitude bigger than I would prefer.

Is there a way to query e.g. OSM for lower-resolution data about borders for all countries?
Alternatively: is there a way to use a tool such as QGIS to simplify borders (hopefully without introducing gaps)?

I'm more than willing to put in the work if needed, but I'm not sure where I should start or which direction I should try. :)


texts: Select specific rows in attributes table by using ArcPy?

texts:
With the zonal statistics as table I made a "park_dem.dbf" file which shows in 28 rows the name of  polygon and also several information about the (max-min-range-area) so my question is,  how can I select the four named specific rows (which I know the name) among the all names in the table and then find that which of these four area has the max range? So after that I have to export this selected area to the shapefile.

So this is my script which I have done. Is there anyone who can help me to solve this?

(I am not a beginner in Python)

from arcpy import env
from arcpy.sa import *

env.workspace = 'ArcCR500_v32/'
arcpy.env.overwriteOutput = True
try:
    if arcpy.CheckExtension('Spatial')=='Available':
        arcpy.CheckOutExtension('Spatial')
    else:
        raise LicenseError
except:
    arcpy.AddMessage('Spatial Analyst License is not available')

#make mix table
in_zone_data = 'ArcCR500_v32/ChranenaUzemi.shp'
zone_field = 'NAZEV_ASCI'
in_value_raster = 'ArcCR500_v32/dem'
out_table = 'Park_dem.dbf'

result_table = ZonalStatisticsAsTable(in_zone_data,zone_field,in_value_raster,out_table, 'DATA')


texts: Pros and Cons on loading vector data from ESRI Shapefile and GeoServer using OpenLayers 5

texts:
I am developing a web map application to show the vector data on top of OpenStreetMap using the open source technologies such as OpenLayers v5.2.0, but am a bit confused to take a decision in what method to choose to render the vector data on map either from GeoServer or directly from ESRI Shapefile. My application will have simple functionalities like feature selection/de-selection,attribute popup etc.

Hence please suggest me the Pros and Cons of these two options to move forward.


texts: How to convert a vector boundary to a raster map?

texts: I need to output a PNG file of a raster map overlayed with boundaries present in a vector map. The problem is, the r.out.png script only exports a raster map. My boundaries are in a vector format. When I tried to convert the vector boundaries into a raster, in order for me to combine them, the converted raster from vector contained the filled vector boundaries. I only need the lines to be present in the raster map. Is there a way to do it? I cannot use the d.rast and d.vect since I am calling the scripts in a shell script so there would not be any GRASS windows open.


texts: How to merge several shapefiles with the same geometries and same attribute names?

texts:
I have data from several runs of a numerical modeling tool. The data always has the same structure - e.g. fixed set of rectangular vector tiles - no matter if there is a "valid" output or zero.
It looks like this (I only put some random data inside):

The attributes also always have the same names - "cell", "value" and "maximum" where "cell" is ID of particular vector tile:

I am probably going to have several tens or hundreds of such datasets and I was asked if it is possible to join them into one vector layer with several attribute columns for each model run.
The ways to merge shapefiles do not expect same geometries so I would probably end with unnecessary overlapping tiles... while the vector tiles remain the same and only the attribute values change.
I would probably be able to put together some CLI linux script first converting all the DBF files to CSV and somehow rename the column names based on the filename to avoid conflicts and finally join it with the geometries - in QGIS, maybe or some sort of geodatabase...
Or is there any better way to do this?
PS: can share some testing data if it helps


texts: Can I use a GeoServer or MapServer service to perform the necessary spatial queries?

texts:
I am building an application and I want it to have the following functionality.

I have a table of points and I would like to be able to query that table for a specific point and get all of the polygons that it is INSIDE of.

For example. If I have a point in Austin,TX, I would like my application to return a reference to a Austin, Texas, United States, North America, essentially grabbing all "parent" polygons. At the same time I would like to be able to do the reverse and say "give me all of the points in the United States", or "give me the points in Texas".

I know how to do the queries with PostGIS, and up until now, I have been loading various shape files into my database from gadm. I have all admin levels, but I am missing American cities so now I have to go download that data. My hunch is I am missing other data besides American cities so I started looking at GeoServer and MapServer.

My question is, am I going about this the wrong way?

Can I query a GeoServer service for the same type of information (giving me all points in a polygon)? Am I wasting my time loading these shapefiles into PostGIS? I am so new to this stuff so any information helps.

PS, right now I am using QGIS to explore the data. Not sure if that helps answer the questions or not


texts:
I'm trying to compile a cython plugin for QGIS.
The plugin is working on windows (32/64) machines but was never used on Mac before.

I compiled the cython for python as described here but there are still issues and code now returns above error.

File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 306, in startPlugin
                plugins[packageName] = package.classFactory(iface)
            AttributeError: 'module' object has no attribute 'classFactory'


texts: Creating list of strings as parameter to Python script tool?

texts:
I am trying to have a list of strings to pass as parameter to my script:

when I start the script it displays:

With "Add Value' button disabled. What am I missing here?


texts: ggmap, ggimage… - Saving and loading maps in R

texts:
In my current R project I'm using the ggmap package to download and display maps. What bothers me is that you need internet access to load maps from google etc. I'd like to load a map from the internet and save it so I can work offline as well. An ideal solution should minimize the amount of code between get_map() (or a related command) and using the saved map within ggmap() (or a related command).

I know already that get_map() has a filename attribute to save maps as an image, and I also know the ggimage() method (I've never actually worked with it, though.) Is there a simple way save such a map, or are these two the best tools I have, and if so, how do you use them effectively?


texts: Adding colours from .style file to multiple shapefiles via script?

texts:
I have added 70 fill colours in a .style file. For adding these colours to my shapefile i would have to manually go to properties/symbology/match to symbols in a style/ then browse to the location where my style file is saved, add it, click on Apply for the colours to change. I do not want to do this for each shapefile.

Is there a script which when run would call this .style file and add its colours to the shapefiles in my workspace?


texts: Multipatch to collada says Background processing has encountered a serious error and is unable to continue?

texts:
I was trying to convert a 2D shapefile of building footprint to 3D SketchUp-readable file.

I first opened the shapefile in ArcGlobe, and in properties set building height as extrusion parameter. Then I used "3D to feature" to convert shapefiles to multipatch.

The last step was supposed to be "Multipatch to collada", one conversion tool, but ArcGIS crashed, saying:

  Background processing has encountered a serious error and is unable to continue

I tried disabling background processing, and then the application crashed when re-running the tool.

Any thoughts what happened, or any suggestions how to convert shapefile (with height value) to collada?

I am using ArcGIS 10.3 student edition, on Windows operation system.


texts: ESRI JS API Can I attach other events to to the map object

texts:
I asked this over on the ESRI forum but it didn't get any traction. I have an app that I'm designing specifically for mobile. At one point in its usage, a user needs to pan and zoom to a location and select a point and get back the evt.mapPoint.x/y. While developing on a desktop/laptop, the double-click event is available. One mobile, it's not. ESRI does support a double-tap/click to zoom the map, but this event doesn't look like I can use it to select a location. I've tried using Dojo Gestures for a tap.hold event, but the ESRI API doesn't honor this event. Right now, this is how I'm implementing it.

pickLocation: function() {
          menus.closePointPicker();
          menus.showTapAndHoldDialog();
          map.disableDoubleClickZoom();
          var pick = on(map, 'dbl-click', function(evt) {
            Event.stop(evt);
            var latLng = webMercatorUtils.xyToLngLat(evt.mapPoint.x, evt.mapPoint.y);
            var location = {};
            coords = {
              latitude: latLng[1],
              longitude: latLng[0],
              accuracy:null,
              heading:null
            };
            var locationInfo = {
              coords:coords,
              point: evt.mapPoint
            };
            pick.remove();
            topic.publish('location-picked', locationInfo);
            menus.hideTapAndHoldDialog();
          });

Any ideas on how to either get a double-tap or tap and hold event to connect?


texts: Calculating Setbacks of Parcel using ArcGIS for Desktop?

texts:
I am looking for any guidance on how to calculate setbacks for a few parcels.

I am using ArcGIS 10.3

I have a shapefile of a few land parcel. They have columns of Side_Setback, Front_Setback, and Rear_Setback.

The image below explains setbacks well.

I am looking to take the original parcel layer and have a model builder and/or python script cut those polygons into a polygon of that parcels buildable area and have the same attributes of the original polygon.


texts: Greek letters of attributes are not shown correctly in FME

texts:
I am using FME 2014 in order to manipulate some shape files.
When I load and inspect the content of the shapes file, I get the following problem:

Although all the greek characters of the content is shown correctly (in UTF-8), the titles of each column (seeing in data inspector) are not recognized.
A similar problem I had some time ago with an excel file which had content in greek. In that case also I could see the content of the excel cells but I could see the names of the tabs (sheets).
Back then it was enough to rename the sheets to latin characters but now I don't have this luxury (some hundreds of shape files).

What is the issue here? Is there some internal encoding properties in FME? I see I have the same problem in FME 2015.


texts: Count of raster cells as variable in ArcGIS ModelBuilder?

texts:
I want to use a cell count of an integer raster as a variable in Arcgis 10 ModelBuilder (and later in a python script), but can't seem to find how to access this information.

Is there a (python) function that will return a cell count for a particular value?
eg. something like

count([raster], value)?

Or how do I access the count column for a particular value in the raster table?


texts: Create polyline layer from csv-file (id, lat, long, time) in QGIS

texts:
Some units are moving and regularly registering it location. Each of these register events corresponds to row in a csv-file. This row contains fields:

id (unit identifier),
lat, long (event location),
time (event time in yyyy-MM-dd HH:mm:ss format, it can be a number of seconds from 1970 if it really matters)

I need to create a polyline layer in QGIS from this csv-file, that represent routes of the units.


texts: System Requirements for 3d GIS & 3d Modelling

texts:
I'm looking for system (hardwares) required for view, edit & analyse 3d GIS & modelling data.
And what is the decent monitor size should be ? single/Dual screens ?

Modern systems comes with 8GB Ram, 2GB Graphic Card, 1TB space, 3 button mouse.
But do we need anything else for a workstation ? which has decent capability for 3d GIS & 3d modelling data ?

And do you reckon 2 different set of systems for viewing and editing capable workstation ?


texts: How to split a single polygon with another polygon (from a different layer) with Qgis

texts:
I have two polygon layers: Zones and Neighbourhoods. I need to manually split one polygon of my Zones layer into two polygons, and the resulting two polygons need to be related to the neighbourhoods.

In the images below, Zones are pink outlined; Neighbourhoods are blue areas. First image is what I have now, second image is what I would like to have (two polygons splitted by the blue line)

I used to do this operation in Mapinfo, following these steps:

making Zones layer editable,
selecting the polygon of the Zones layer to be cut,
then Objects > Set Target,
then select the polygon of the Neighbourhoods layer that serves as the clipper,
Object > Split,
Clear Target.

Can I do something simmilar in Qgis (I am using version 2.18.11)?


texts:
I'm using esri-leaflet, and leaflet inherently to add a point feature class (geojson) to a map. What is the best way to add and style the layer? I've see pointToLayer on the leafletjs.com, but don't want to use the circleMarker. Below is a code snippet of what I am working with at the moment:

new L.esri.featureLayer(schools, {
    style: function (feature) {
        return{
            weight: 2,
            opacity: 1,
            'marker-color': '#FFF'
        }
    },
    onEachFeature: function (feature, layer) {
        layer.bindPopup(L.Util.template(m_template, feature.properties));
    }
}).addTo(map);

I would like to use AwesomeMarkers, but I don't know how to implement that with a featureLayer. From my understanding featureLayer extends L.GeoJson. The 'marker-color' has not worked. A live preview of this is available HERE


texts: I have just downloaded Qgis2threejs plugin from: Plugins/Manage and install plugins but the tool doesn't appear in the toolbar. Does anyone know what is wrong? I'm using QGIS 2.14.1 Essen.


texts: Calculate distance to points along line in ArcMap?

texts:
I have two shapefiles: a series of points and a polyline created from the points with the tool Points to Line. Both shapefiles have the same projected coordinate system with a linear unit in metres.

If you were to calculate the distance to the points along this line - which is what I'm trying to do - the first point would have a value of 0 metres and the last a value of closer to 40 000 metres (the length of the line). I want to gather these values in an attribute field or a table so that I can plot the points in a graph with the distance in metres on the x axis.

I have tried the tool Locate Features Along Routes but for some reason the data doesn't come out right. I have tried the tool with both the original line and a route created from the line with the tool Create Routes.

Input values (identical for line and route):

Results - line:

Results - route:

Do you have any idea where I'm going wrong and what I should be doing instead? Is there another way to go about this?

I use ArcMap 10.4.1.


texts: In ArcGIS for Desktop there is an option to graphic clip inside the data frame - this clips the feature you're viewing without actually clipping the feature into a new feature.  Does Quantum GIS have a tool to do this?


texts: Importing EXIF Metadata into attributes fields in QGIS

texts:
I use to imports JPG images with "ImportPhotos". It lacks some importants EXIF Tags such : ISO , Shutter Speed, aperture, resolution and Field of view.

There is an answer about how to import XMP Metadata using the layer created by ImportPhoto, using the path (stored in attributes fields).

I need help to complete this script to be able to import as well the EXIF TAG/values listed above.

You can find some sample images there: https://we.tl/t-BjPuzBGMRc