texts: Publishing Google Sheet to web as GeoJSON file?

texts:
I am trying to create a map in Mapbox to be embedded on a website and have the embedded map automatically update every time the source data is updated. I have my source data in a Google sheet and have discovered how to manually convert that into a GeoJSON file and then manually import that file into Mapbox, so that the data populates my map. However, I would like to find a way to automate this process. My proposed solution is that I publish the sheet to the web as a GeoJSON file and then set it to automatically republish whenever changes are made. However, I don't see an option for publishing a sheet as a geoJSON file.

Does anyone know of a script for this?


texts: I use ogr2ogr to import an .sqlite file into my PostgreSQL database. The new database table I get has a column of type "character varying". However I want it be of type integer as all the values in the particular column are integer numbers with the exception of some empty fields. How do I do this?


texts: Looking for a way to change the default shapefile label in ArcGIS 10

texts:
Is there a way to change ArcGIS's default shapefile label Font: Aerial Size: 8 to something we are using over and over again? As shown in attached screen shot, I'd like to include a label expression in addition to a saved text symbol:

It would be very helpful to save the couple of steps it takes to change labels as I am working with hundreds of maps / shapefiles.


texts: How to limit/page through results from ArcGIS SearchCursor?

texts:
I would like to get for example, the first 200 rows from a database by ArcGIS Python SearchCursor. How can I do that? I've got a featureclass which has about 500k features and reading it all in one shot is taking to long. My present code is:

fields = ";".join([field, "OBJECTID"])

rows = gp.searchcursor(connectionString, "", "", fields, "OBJECTID")
row = rows.Next()
while row:
    values.add(row.getValue(field))
    row = rows.Next()

return values

I've tried to set counter and stop iterating when it came to 50, for example. But when I would like to take another 50 records I don't know how to start from the last one that i fetched.

fields = ";".join([field, "OBJECTID"])
counter = 0

rows = gp.searchcursor(connectionString, "", "", fields, "OBJECTID")
row = rows.Next()
while row:
    counter += 1
    values.add(row.getValue(field))
    row = rows.Next()
        if (counter >= 50):
            break;

return values


texts:
I read How can I label the coordinate position of the grid coordinate lines but I still do not understand how to use the LineFeatureBuilder class to label the grid lines.  Could someone let me know how do I use this class in my code to get the labels?
Code creating the lines:

        List<OrthoLineDef> lineDefs =
                Arrays.asList(
                        new OrthoLineDef(LineOrientation.VERTICAL, 2, 10.0),
                        new OrthoLineDef(LineOrientation.HORIZONTAL, 2, 10.0)
                );
        double vertexSpacing = 0.1;
        SimpleFeatureSource grid = Lines.createOrthoLines(gridBounds, lineDefs, vertexSpacing);
....
        Style latlonStyle = SLD.createLineStyle( strokeColor, strokeWidth );


texts: How to export Atlas as images from a headless server?

texts:
I have a QGIS project that has an Atlas print composer that I can use to generate multiple images (JPG) of each of my departments.

Now I would want to generate all those images but directly from a Linux server which does not have any graphical server (it is a headless server, like many Linux servers).

Is this possible?


texts: Can you assign and use variables within text box QGIS print composer

texts:
I'm using dynamic text within a text box in print composer in QGIS 3.12.  Is it possible to assign and use variables within this to avoid repeating the same call several times?
I'm using [% attribute(@atlas_pagename, "Acres) %] to pull up the area of the current atlas feature, then later using the same number to calculate the percent of different features.  Is it possible to assign the value of that call to a variable and use it later?


texts: Can GPS transmission time be the same for two different ephemerides?

texts:
I am seeing a RINEX navigation file where I have two different ephemerides for the same satellite at different times. Nothing surprising. The thing is that they both have the same transmission time, which I understand it is derived from the Z-count from the Handover word. Wouldn't the transmission time need to be different?

I'm not totally familiar with the GPS navigation message. Is it possible that a HOW is sent, then a word containing a GPS ephemeris for a particular SV and then a new word containing ephemeris (different IODE) for the same SV sent again? Is it GPS ICD compliant?


texts:
I am creating trying to create a function inside PostgreSQL from python. My connection to PostgreSQL database uses psycopg2 and connects successfully in other instances. Code:

pg_cursor.execute('CREATE OR REPLACE FUNCTION %s.fix_geometry(geometry)\
                  RETURNS SETOF geometry AS\
                $BODY$\
                        SELECT geom the_geom\
                        FROM \
                                (\
                                SELECT (st_dump(st_buffer(st_snaptogrid(st_makevalid(the_geom), 0.5), 0))).geom\
                                FROM (SELECT geom the_geom FROM st_dump(st_snaptogrid($1, 0.5))) a\
                                ) b\
                        WHERE geometrytype(geom) = \'POLYGON\' AND st_area(geom) >= 0.01;\
                $BODY$\
                  LANGUAGE sql VOLATILE\
                  COST 100\
                  ROWS 1000;' %(schema)) #line 84
pg_connection.commit()
pg_cursor.execute('ALTER FUNCTION %s.fix_geometry(geometry) OWNER TO analysis' % (schema))
pg_connection.commit()

I get error:

line 84, in
    ROWS 1000;' %(schema))
ProgrammingError: type geometry does not exist

When I run the code in pgAdmin it executes successfully. What am I missing?

I had originally put this question here: https://stackoverflow.com/q/41602531/4641482 but have yet to have had a response, maybe wrong forum?

Python 2.7, PostgreSQL 9.3


texts: Installing Spatialite on Centos7: select mod_extension() doesn't work

texts:
I'm trying to install Spatialite on a CentOS 7 server.

Sqlite (version 3.7.17 ...), works fine.

Using the command

yum list | grep spatial

the result is the follow

So I've installed Spatialite using the command

sudo yum install libspatialite

When I try to check abot spatialite version I obtain ...

If I search for "spatial" libraries in my system I obtain ...

Any suggestions?


texts: relative path for image elements in ArcMap

texts:
Is there a way to set relative path for image elements (like a company logo) in ArcMap projects?
The classic "Store relative path" option in the Map Properties doesn't seem to be working. Thanks!


texts: Reducing a MapInfo region to within a radius

texts:
I'm trying to get rid of parts of a region outside of a 20km radius. At this point I have a region A and a circular region C that are centered at the same point.

I tried using overlap(A.obj, C.obj) and it kept crashing, so I tried using check regions to check for overlaps, but it returns the statement:

  "Check Region did not find any data problems"


texts: Transfer road netwrok node id's from shapefile to .osm file using ogr2osm

texts: I have been able to convert a custom road network in shapefile format to .osm data. The original shapefile contains a column "source" and a column "target" to indicate the node id's. After the conversion to .osm format however, these node id's appear to have lost. Instead, new node id's have been chronologically numbered in the .osm file (see below). Is there a way to transfer the original node id's from the shapefile to the .osm file?


texts: How to measure length of line accurately?

texts:
In postgis 1.5 for the below given query the answer is 1540651.62674101 meters ,

SELECT ST_Length(ST_GeomFromText('LINESTRING(77.83042478941326
9.98463896136159,78.03469506809812 23.90398909552571 )',4326),true);

but in google earth the same line gives 1538923 meters in its properties, around 0.1% error between the two readings.

Which one is accurate, if Google earth reading is accurate, how to get it in postgis1.5 ...

any idea...

please help.


texts:
I'm trying to compile MapServer 7 with:

Proj 6
Geos 3.7
Gdal 2.4
PostGIS 2.5
PostgreSQL 10

But for some reason it complaining about pg_config not found, I usually insert the path /usr/pgsql-10/bin/pg_config to compile software that need the same file but with MapServer 7 I can't

Here my cmake :

cmake -DCMAKE_INSTALL_PREFIX=/opt \
        -DCMAKE_PREFIX_PATH=/usr/pgsql-10/bin:/usr/local/src/postgis-2.5.2:/usr/local:/opt:/usr/local/src \
        -DWITH_CLIENT_WFS=ON \
        -DWITH_CLIENT_WMS=ON \
        -DWITH_CURL=ON \
        -DWITH_SOS=ON \
        -DWITH_FCGI=0 \
        -DWITH_PERL=0 \
        -DWITH_RUBY=0 \
        -DWITH_JAVA=0 \
        -DWITH_CSHARP=0 \
        -DWITH_PYTHON=0 \
        -DWITH_SVGCAIRO=0 \
        -DWITH_FRIBIDI=0 \
        -DWITH_PROTOBUFC=0 \
        -DWITH_CAIRO=0 \
        -DWITH_HARFBUZZ=0 \
        -DWITH_ORACLESPATIAL=0 \
        -DWITH_MSSQL2008=0 \
        ../ >../configure.out.txt

I don't really understand why I can compile other software correctly and MapServer no, I don't even understand what should I put in the path as looking at this code (supposedly responsible for the check)

https://github.com/Kitware/CMake/blob/master/Modules/FindPostgreSQL.cmake

this is the first error that appear:

Cmake Warning at cmake/FindPostgreSQL.cmake:21 (message)
Pg_config not found, will try some defaults
Call stacks (most recent call first):
cMakeLists.txt:755 (find_package)

and the second one is:

cMake error at CmakeLists.txt:73 (messages)
POSTGIS library/components/dependency could not be found.

HINTS:

-disable postgis...
- add cmake_prefix_path...
Call stack (most recent call firsts)
CMakeLists.txt:590 (report_optional_not_found)

How can I compile it with PostGIS support?


texts: Correct planet radius, for calculating metres per pixel?

texts:
I am putting together some code to create tiles from various maps, for various planetary bodies, but I am not sure which value of planet radius I should be using, for indicating the metres/pixel during display?

When I look at the Open Street Maps entry on 'zoom levels' a value of 6372.7982 km is given.
Looking at the 'Earth Fact sheet' I see:

Equatorial radius (km):         6378.137
Polar radius (km):               6356.752
Volumetric mean radius (km):     6371.008

Wikipedia provides a value of 6371 km, which corresponds to the 'Volumetric mean radius' of the NASA fact sheet.
Google appears to be using the value of 6378.13 km, basing on the circumference indicated in this other Stack Exchange answer.

Ideally I would like to use the data from the 'Planetary Factsheet' (for consistency across planetary bodies), so I am wondering whether I better off using the 'Equatorial radius' or the 'Volumetric mean radius' as reference value? BTW I should note that I will be using this along with metres per pixel equation, specified on the Open Street Map entry, on 'zoom levels':

S=C*cos(y)/2^(z+8)


texts: PostGIS is rejecting an srid code for my projection. I've found a nearly identical projection w/ a legit srid. Will it work?

texts:
My projection is this: NAD 1983 StatePlane North Carolina FIPS 3200 Feet, which has a proj4 string that looks like this:

+proj=lcc +lat_1=34.33333333333334 +lat_2=36.16666666666666 +lat_0=33.75 +lon_0=-79 +x_0=609601.2199999999 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs

Its SRID is 102719. I've tried to create a db column in PostGIS with an srid option for this projection, but its rejecting the srid, calling it invalid.

I've found a similar projection with a valid srid code, NAD83 / North Carolina (ftUS) (srid=2264), which has a proj4 string nearly identical to the above:

+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.2192024384 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs

As you can see, the only difference between these two is that the lat_1 and lat_2 degrees are swapped. They even use the same three degrees, only the first two are differently ordered.

Can I use this second projection?

UPDATE: The answer seems to be leaning towards 'yes', I just ran projection-to-projection transformations from a single projection to the NC projection using the "legitimate" and "illegitimate" proj4's and got similar results down to the 7th decimal place


texts: Why aren't my CartoDB map layers loading?

texts:
We have built our own instance of CartoDB on our servers, setup production environtment successfully. We can upload shape files, table and map view works great.

When we trying to create visualization with two different data and publish, CartoDB generates a short url. But when we try to access the visualization (sharing) through that short url none of the layer gets displayed on map. The basemap loads fine, but the layers don't appear.

I could see the following error message on console.

Failed to load resource: the server responded with a status of 404 (Not Found)  http://cartodb.com/noneuser.html

I really don't know what's wrong, could someone can help me on this? Thanks in advance.


texts: How to display lines with Arrows on a FeatureLayer?

texts:
I've setup a feature service with a symbology.  One symbol is a line with an arrow at the end.  The description for such a symbol is not part of the FeatureServer response data.  So it is not supported?  I only see this (not sure why the width is 6):

UniqueValueInfos:
Value: Arrow
Label: Arrow
Description:
Symbol:
Style: esriSLSSolid
Color: [0, 0, 0, 255]
Width: 6

If not supported, any tips on how I might proceed?  Perhaps a termination icon attribute and then some custom gfx code?  Anything already exist for this?


texts: Using Polygon centroids tool - points aren't always inside polygon, can I fix that?

texts: When I use Polygon centroid tool, points aren't always inside polygon, because some polygons have weird shape with holes inside them. I have big amount of polygons and I need to use this tool and I need all centroids to be inside polygons. Is there any way how can I make them to be inside polygons?


texts:
Is there a way to locate diagrams and move them wherever you want in QGIS 1.7.4/ 1.8.0?
I created a pie chart under the tab "diagram". I can mark those diagrams with a green box and move them via the label-tool, but they won't change their position in the end. When I create diagrams with the "diagram overlay" tab they cannot even be marked.

Also, what is the difference between "diagram" and "diagram overlay"?


texts: How to select the polygons based on a raster?

texts:
I have many rasters and one vector layer with many polygon features and I am looking for a way to find the polygons which fall above each raster (sort of select by location)

Please suggest an way so that I can automate the process.

End Goal is to find number of polygons intersecting each raster and their names.


texts: 3D polylines, recieving z-value as attribute in ArcGIS for Desktop?

texts:
I am using ArcGIS for Desktop, and attempting to create a runoff mode for a water course and the surrounding floodplain.

I have depth data in the form of contours, given as 3D-polylines. Each line has an attibute which shows the depth at each line but i need the z-value instead since i need to merge the depth contours with an existing elevation model outside of the water course. In order to do that i have to have the elevation as an attribute (a z-field in the attibute table) for each line. How can give each line the Z-value as an attribute?

When i extract the vertices i use add xy-data which gives me the z-value. I wish to do something similar but for the polylines. Each lines has only one z-value assigned to it which is clear when reviewing the z-value in each vertex.


texts: I want to import a table (created and formatted in MS word?) to appear as text in the MapInfo layout window.  Obviously this isn't individual data records which appear in the MapInfo browser window.  So would I need to format the MS Word table, save the file as a JPG/png format then import into MapInfo as a raster file?  Is there another way in MapInfo of doing this - just type the table contents as text into the layout window then draw a rectangle around the text and make this box transparent but with a border?


texts: How to split rasters tile size not in pixels but split in meters?

texts:
I write my raster calculator, similar in functionality to Qgis, only I will read arithmetic operations with Pyopencl
I'm stay in only one place, there are two raster images that have the same size in meters, but a different pixel resolution
raster1 =0.037368,-0.037368
raster2=0.0463053,-0.0463053

Since my rasters occupy more than 1 GB ,I'm forced to break them into buffers 512 to 512

bufferXSize = 512
bufferYSize = 512
for x in range(numX):
    col = 0
    print col
    for y in range(numY):
        xsize = rasterXSize - x*bufferXSize if x == numX-1 else bufferXSize
        ysize = rasterYSize - y*bufferYSize if y == numY-1 else bufferYSize
        slice = gdalData.ReadAsArray(x*bufferXSize, y*bufferYSize, xsize, ysize)
        r1 = [rx, ry, rx + (geoT[1] * slice.shape[0]), ry + (geoT[5] * slice.shape[1])]
        rx+=slice.shape[0]
        ry+=slice.shape[1]

My question is, because the two images have the same size in meters, how can I not retreat by 512 to 512 pixels, and what would the indentation of blocks be done in meters (for example, indent 25 to 25 meters on the first raster, and the same indent on the second)


texts: DIfference between Thumbnail and GeoTIFF Images in GEE

texts:
I am trying to build a Deep Learning model that takes Satellite Image as Input and uses a few demographics parameter as output.

I am new to Google Earth Engine, so I was checking different ways to export the data.
I can export it as a GeoTIFF Image or a Thumbnail. Since I just need the RGB bands of the Image, will it make any difference if I use a thumbnail.
For eg: Will I lose any information that might be helpful for the model to learn?


texts: Trying to extract point X and Y from Geodatabase from outside a GIS

texts:
I have a geodatabase which I would like to integrate with processes outside ArcGIS.  (On our systems Arc is a fairly closed system).

I have managed to connect to the database and perform an extract using OLEDB but cannot get the X,Y coordinates directly. Instead there is a SHAPE field which contains binary data:

SHAPE
346C0000010C80EB51B88BF81B4152B81E8539A11041
346C0000010CD0A3703D79B41B418FC2F5287A381041
346C0000010CF051B81EB7651B4166666666E50C1841
346C0000010C6066666655061C4148E17A14A4A51041
346C0000010C10AE47E176C31D4185EB51B8423B1641
346C0000010C703D0AD7D6BB134166666666C0C41241
346C0000010C5050FC98D4B813416FF085C96BC61241
346C0000010C400AD7A31C3615417B14AE4790D91141
346C0000010C2085EB5139591441B81E85EB1BAE1341
346C0000010C703D0AD7BA671541CDCCCCCCD50D1441
346C0000010C608FC2F5A36015410AD7A37031881441
346C0000010CF051B81E25211441B81E85EBC62B1441
346C0000010C2085EB51E382154166666666F8211441
346C0000010CC0F5285C1745184152B81E8505B91041
346C0000010C2085EB516B961741713D0AD7CBE91441

Does anyone know how this binary data should be interpreted? I assume it is a struct, 1 part declaring it is a point, and 2 other parts declaring the X and Y?


texts: Add fields from a feature to a feature in another layer

texts:
I have a criteria to select the features from a layer that overlap the features of another layer, but i do not know how to copy the fields at the same time:
Tyhe code does:
- two layers (1 and 2) of polygons
- if one feature of polygon2 overlap more than 80% one of polygon1, it selects it.

How to also copy the fields of polygon1 to add them to the fields of the selected polygon2?

from qgis.core import *
import processing
import os

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

#polygon1
feats_lyr1 = [ feat for feat in layers[0].getFeatures() ]

#polygon2
feats_lyr2 = [ feat for feat in layers[1].getFeatures() ]

selected_feats = []

for i, feat1 in enumerate(feats_lyr1):
    for j, feat2 in enumerate(feats_lyr2):
        if feat1.geometry().intersects(feat2.geometry()):
            area1 = feat1.geometry().intersection(feat2.geometry()).area()
            area2 = feat2.geometry().area()
            crit =area1/area2
            if crit > 0.8:
                selected_feats.append(feat2)

epsg = layers[0].crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&index=yes"

mem_layer = QgsVectorLayer(uri,
                           "mem_layer",
                           "memory")

for i, feat in enumerate(selected_feats):
    feat.setAttributes(selected_feats[i].attributes())

prov = mem_layer.dataProvider()
attr = layers[1].dataProvider().fields().toList()
prov.addAttributes(attr)
mem_layer.updateFields()
prov.addFeatures(selected_feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)


texts: How to copy/move Mosaic Dataset to another location?

texts:
On my previous question on Mosaic Dataset vs. ERDAS 2011 Is ERDAS 2011 or ArcGIS 10.1 best to use when mosaicking 100's of ECW image files?, I did research for a bit was wondering what happens if you use the Mosiac Dataset and decided to copy and move to another computer. It would display a black and white checkerboard.

Do I need to copy the original ECW files along with the Mosaic Dataset geodatabase for it to make it work or display it ?

When I load the moasic dataset to ArcMap, the only way to see the imagery is to zoom in ?
Is there a specific scale that can I see one set of imagery ?

Any suggestion ?

One problem with ERDAS is that I can do it but the image is actually came from ArcSDE we recieved from BLM and when I review it and it has the tiles which doesn't show overlaid each other. This is the reason why I choose Mosaic Dataet.


texts: R - Occasional error when using alphahull package areaahull() function

texts:
I've recently started using the fantastic alphahull R package, which works beautifully most of the time but occasionally I am getting the following error when running the areaahull() function:

Error in row:rownew : NA/NaN argument

For most of my data the function works as expected, but with some data and for certain values of alpha, the above fatal error occurs. It does not seem to depend on the complexity of the data, sometimes the error occurs with a few dozen data points, while it will run fine with a dataset of a few thousand points.

The code is very straightforward, simply by changing the value of alpha for a given dataset I can generate the error.

library(alphahull)

lonlat <- read.csv("areaahull_error.csv")

alpha <- 0.2

alphahull <- ahull(lonlat, alpha = alpha)
areaahull(alphahull)

I'd be interested to know if anyone else has seen this error and if there's anything to be done about it, other than trying to catch it, or tweak alpha to avoid it when it occurs.


texts: Why is my field not nullable when no value is NULL?

texts:
I am quite new to arcpy.
My script seemed to work just fine but I get an error, "the field is not nullable". I guess the problem is with the last block of code but I don't understand why nor how to get around.

Some of the original columns ['T2TC1', 'T2TC2', 'T2TC3','T2TC4','T2TC5'] are equal to 0. No value in my "TOT_Count" column is equal to 0, so I am not dividing by 0.

TOT_Count is the sum of ['T2TC1', 'T2TC2', 'T2TC3','T2TC4','T2TC5']. This seems to be creating the problem but I don't know why. If I replace the values of this column by a random number the problem disapears.

import arcpy
table = "C:\Users\Me\Documents\datas\Count\FinalT4.shp"
arcpy.env.overwriteOutput = True

#input fields with counts#
fields = ['T2TC1', 'T2TC2', 'T2TC3','T2TC4','T2TC5'] # etc.
newFields = []

#create a new f field corresponding to each field#
for field in fields:
    field_name ="F_" + field
    arcpy.AddField_management (table, field_name, "FLOAT")
    newFields.append(field_name)

#for each new field created, calculate the old field divided by the sum#
for newField in newFields:
    for field in fields :
        expr =  '(!{}! / !TOT_Count!)'.format(field)
        arcpy.CalculateField_management (table, newField, expr, 'PYTHON')


texts: Dissolving roads of same levels in ArcGIS Desktop?

texts:
I have a 3D road shape files. Roads are of different levels, highway, primary,secondary and so on. And when i used the tool feature to line i splited every line in every intersection.

Now i want a way to dissolve roads that are snapped at intersection and have the same level (z value). In this way they do not dissolve with other road that pass on the same intersection but on different level (z value).

Like this picture:


texts: Why does reading shapefile from PostgreSQL database give Error reading OID from table ...?

texts: I'm using arcgis 10.1 and i have a connection to a PostgreSQL 9.3 database. When i try to load a shapefile using the add query layer I get an error when opening the attribute table. The error says "Error reading OID from table. Reading rows has been stopped. Check that the datasource is valid. OID mapped column has null values. the operation is not supported by this implementation". I know there is a field that is causing this issue yet i don't know which one is or how to go solve it. Has anyone come across to the same issue?


texts: Line direction histogram plugin vector identification?

texts:
I was wondering if there is a way of identfying which vector ids fall into a bin. Is this possible to find out? If there isnt is there another plugin that allows you to identify what the angle in a similar fashion as the line histogram. I have used this plugin as a visual aid but do require to know which vectors fall into i.e. 0-11.25 degrees for stats purposes

Below if the line direction histogram and the CSV file that accompanies it..... I want to identify which vectors go in which bin for example the first line on the csv file gives you the sum length and the number of trails in that bin. But I want to identify which vector falls in a bin to carry out statistical analysis


texts: Creating fine interval contour layer from DEM using QGIS?

texts:
I would like to create 20 meter contours for a large region, roughly 120,000 square kilometers.  DEM data is available as TIFF or CSV throughout this region at 0.75 arc second resolution.

What tools should I use to create a coverage of 20 meter contours throughout this region?  I am concerned that large file sizes will present a problem.

I have tried using the contour tool (Raster>Extraction>Contour).  As well as the Contour Plugin which uses DEM point data as an input.  I am running into trouble with large file sizes.


texts: Esri JavaScript API ApplyEdits update not working

texts:
When I attempt to applyEdits as an update to my application I receive the following errors:

I uploaded my most recent revision on Github at: https://github.com/csergent45/streetSigns/tree/119dc2ec277a347c066f92fa13e0ca76f182fb8f

My code is is in app/main.js

You can download the app and run it; don't allow for geolocation and then click on a red point; this will open the sign editing form; you can then modify it, click update and then an error will occur in the console where applyEdits is on line 723 in the app.


texts: How to query all layers in ArcMap and get access to them by name or index?

texts:
I have ArcMap with tons of layers. Some of these layers may contain another layers with data (aka Composite layers) or they don't. The depth of composite layers in unknown. I need to query all layer names and put it to comboBox (done). Next step is finding selected in comboBox layer (by name? by index?) and get all its fields to another comboBox.

My problem is I can't get access to layer by name or index because they may be represented as composite layers and my method for finding layers by name or index will not work. How should I get direct access to layer or maybe create custom array with indexes and names?

Here is my code:

private List<ILayer> GetAllLayers(IMap mp, ICompositeLayer cl)
{
    listIndexesAndNames.Add(new comboBoxLayerListItemsClass(){ItemName = "", ItemIndex = 0} );

    if (mp == null && cl == null)
    {
        return null;
    }
    //ILayer l;
    List<ILayer> listOfLayers = new List<ILayer>();
    if (mp != null && cl == null)
    {
        for (int i = 0; i < mp.LayerCount; i++)
        {
            if (mp.Layer[i] is ICompositeLayer)
            {
                cl = mp.Layer[i] as ICompositeLayer;
                GetAllLayers(null, cl);
                ++index;
                listIndexesAndNames.Add(new comboBoxLayerListItemsClass() { ItemName = mp.Layer[i].Name, ItemIndex = index });
            }
            else
            {
                listOfLayers.Add(mp.Layer[i]);
                comboBoxLayerList_Houses.Items.Add(new comboBoxLayerListItemsClass() { ItemName = mp.Layer[i].Name});
            }
        }
    }
    else if (mp == null && cl != null)
    {
        for (int i = 0; i < cl.Count; i++)
        {
            if (cl.Layer[i] is ICompositeLayer)
            {
                ICompositeLayer cl2 = (ICompositeLayer)cl.Layer[i];
                GetAllLayers(null, cl2);
            }
            else
            {
                listOfLayers.Add(cl.Layer[i]);
                comboBoxLayerList_Houses.Items.Add(new comboBoxLayerListItemsClass() { ItemName = cl.Layer[i].Name });
            }
        }
    }
    return listOfLayers;
}


texts: I am using QGIS 3.8. I am creating a geological map, and because I wanted to merge two different polygons, it created a separate layer. I want to convert this dissolved layer into the attribute to my geology layer with the other polygons. Sorry for my english termology, I'm not an experienced user. Can you help me?


texts:
Below is my code for calculating NDVI using Landsat 8. I have used band 4 and 5, but I am not getting any output.
How should I get the output image?

import numpy as np
import rasterio
import subprocess
b4 = rasterio.open('C:\\Users\\ARSH\\Downloads\\LC08_L1TP_143051_20191219_20191226_01_T1\\LC08_L1TP_143051_20191219_20191226_01_T1_B4.tif')
b5 = rasterio.open('C:\\Users\\ARSH\\Downloads\\LC08_L1TP_143051_20191219_20191226_01_T1\\LC08_L1TP_143051_20191219_20191226_01_T1_B5.tif')
red = b4.read()
nir = b5.read()
np.seterr(divide='ignore', invalid='ignore')
# Calculate ndvi
ndvi = (nir.astype(float)-red.astype(float))/(nir+red)

# Write the NDVI image
meta = b4.meta
meta.update(driver='GTiff')
meta.update(dtype=rasterio.float32)

with rasterio.open('NDVI.tif', 'w', **meta) as dst:
    dst.write(ndvi.astype(rasterio.float32)


texts: Performing complicated calculation on bands in Google Earth Engine

texts:
I'm trying to perform ANIR index calculation on my image collection. I'm having trouble writing the code in the right way. ANIR index formula is :
cos^-1[((a^2)+(b^2)+(c^2))/(2*a*b)]

(a2, b2 and c1, are Euclidean distances between bands B4, B8, B11/B12 respectively).
This is the function I have written so far-
var addANIR = function(image) {
  var ANIR = image.expression(
  '(1/cosin)',
  {
    'B4': image.select('B4'),
    'B8': image.select('B8'),
    'B11': image.select('B11'),
    'A':((ee.Number(832.8-644.6).pow(2)).add(image.select('B8').subtract(image.select('B4')).pow(2))).sqrt(),
    'B':((ee.Number(1613.7-832.8).pow(2)).add(image.select('B11').subtract(image.select('B8')).pow(2))).sqrt(),
    'C':((ee.Number(1613.7-644.6).pow(2)).add(image.select('B11').subtract(image.select('B4')).pow(2))).sqrt(),
    'top':A.pow(2).add(B.pow(2)).add((C).pow(2))/ (ee.Number(2).multiply(A).multiply(B)),
    'cosin': top.cos()

  }).rename('ANIR');
  //NDVI index

  return image.addBands(ANIR.add(ee.Image(0.0)));
  //Add 0.0 value to NDVI (currently not in use)

 };

This function gives me an Error. What is wrong with the function/code?


texts: Creating new mapset with Python and EPSG 4326 (Invalid north field 1N)

texts:
I am using grass 7.8. I have been trying to create a new location and then a new mapset in python. It works fine for any EPSG code except 4326. With EPSG 4326, when I try simple commands (like printing the region for instance) it raises an error:
"ERROR: invalid north field 1N"
I noticed that what was causing the problem is the WIND file which contains the region of the mapset. When the ESPG is 4326 then the north field gets as value "1N". When I create the mapset and location directly with GRASS GIS and without using python I have no problem in this case even though the WIND file looks exactly the same when I create the mapset with python. Does anyone have an idea about this? This is my code:
grassbin='C:\\\"Program Files\"\\\"GRASS GIS 7.8\"\\grass78.bat'
os.environ['GRASSBIN']= grassbin
from grass_session import Session
import grass.script as gscript

import grass.script.setup as gsetup
with Session(gisdb=mygisdb, location=mylocation, create_opts='EPSG:4326'):
     print("Created new location")

with Session(gisdb=mygisdb, location=mylocation, mapset=mymapset, create_opts=''):
     print("Created new mapset")

gsetup.init(gisbase.decode(), mygisdb, mylocation, mymapset)
gscript.run_command('g.region', flags="p")
gsetup.finish()


texts: How to know if a coordinate is within a polygon in shapefile not working?

texts:
I have

coordinates (lat and lng) in an Excel document
a shapefile which contains all the different Canada provinces shapes.

I would like to be able to generate a new field in Excel in order to classify the different coordinates into the different canada provinces, census division and sub division

I tried the below code, but it is not working.

import fiona
import shapely.geometry

with fiona.open(r"D:\Users\Jonathan\Desktop\CRA-Project v2\Census Division\lcd_000b16a_e.shp") as fiona_collection:

shapefile_record = fiona_collection.next()

# Use Shapely to create the polygon
shape = shapely.geometry.asShape(shapefile_record['geometry'])
#print(shape)

point = shapely.geometry.Point(46.362914,-63.503809) # longitude, latitude

# Alternative: if point.within(shape)
if shape.contains(point):
    print("Found shape for point.")

Update 1:

point = shapely.geometry.Point(46.362914,-63.503809)

Polygon:
Link

Update 2:

import fiona
import shapely.geometry

with fiona.open(r"D:\Users\Jonathan\Desktop\CRA-Project v2\Census Division\lcd_000b16a_e.shp") as fiona_collection:

shapefile_record = next(iter(fiona_collection))

shape = shapely.geometry.asShape(shapefile_record['geometry'])
#print(shape)

point = shapely.geometry.Point(573339.364,6177352.077)

if shape.intersects(point):
    print("Found shape for point.")

UTM:

point = shapely.geometry.Point(573339.364,6177352.077)

New Polygon: Here


texts: CartoDB - random points being added to map

texts:
I'm working on a tool that allows users to draw geometries on a Leaflet map, which are then sent to CartoDB.  However, I've noticed an odd behavior... after a polygon has been added, a varying number of seemingly random points get added to the map as well.

For example, I drew a rough polygon of Manhattan, but here's what the map displays...

What strikes me as particularly odd is that, if I remove the polygon CartoCSS, the points remain... some of them can be clicked on and the little info window that pops up correspond to the ID of the polygon, but others are unclickable.

For reference, here is the data from the_geom column for the Manhattan polygon...

{"type":"Polygon",
 "coordinates":
               [
                [
                 [-74.00802612,40.75037809]
                 [-74.01626587,40.70250471],
                 [-73.9956665,40.70979201],
                 [-73.9805603,40.71499674],
                 [-73.97232056,40.74517613],
                 [-73.93112183,40.79821703],
                 [-73.9352417,40.8335541],
                 [-73.91052246,40.87302619],
                 [-73.92837524,40.87614141],
                 [-74.00802612,40.75037809]
                ]
               ]
}

And I checked these... not one of them corresponds to these mystery points.

Does anyone know why this is happening and what I can do about it?


texts: Using ArcGIS in mixed version environment?

texts:
What are some best practices for ArcGIS 9.3 and 10.0 in the same shop?

We have a couple of dozen users on 9.3. The bulk of our data is in file geodatabases and geotiffs, with a respectable splash of shapefiles and occasional jaunts to ArcSDE (9.3.1) hosted by another department. We will be introducing 10.0 to the power users first and for an unknown period of time, likely 6 months, the two will be operative at the same time. All of the data as well as map compositions, layer files, symbol styles, etc. will be used by both groups.

How do you structure things so these two sets can peacefully co-exist? What are the things to watch for?


texts: ESRI Silverlight Template - 'View Designer' of MainPage.xaml Crashes Visual Studio

texts:
At last, first adventure with the ESRI Silverlight API...

So, I have installed ArcGIS API for Silverlight v2.4 and all its dependencies.  Then, I downloaded the 'Showcase' Template.  Opened the project in Visual Studio Ultimate 2010, and the project built and compiled fine (with some warnings about ScaleBarUnit).  Hit F5, and application ran fine out of dev server.

Then, went to go look at xaml markup in MainPage.xaml by right clicking and choosing 'View Designer' option.  Xaml starts loading, then Visual Studio crashes.

I've spent the last hour or so googling the issue, and haven't found the fix.  I am hoping someone may know what  I need to do to view the xaml markup in Visual Studio designer?


texts: heat map plugin from QGIS python console

texts:
I am trying to use the heat map plugin from the processing console. On the QGIS documentation it says the plugin uses Kernel Density Estimation. I have used the saga Kernel Density Estimation but do not get results that look similar to the heat map.
 How can I import the plugin into python console.
 For vector it's from fTools import *.
How can I use  raster plugins in processing i.e heat map plugin and georeferencer.


texts: Giving user option to show report on screen or create text file using Python Script Tool of ArcPy

texts:
I have been working on this script tool in ArcGIS Pro. The following code is what I got so far to create a tool in ArcGIS Pro, connected to the script, that takes parameters from users (initially, the geodatabase filename) and shows the report that lists the feature classes it contains to the console  on the screen.

import arcpy
import os
arcpy.env.workspace = arcpy.GetParameterAsText(0)
featureclasses = arcpy.ListFeatureClasses()
for fc in featureclasses:
    arcpy.AddMessage(fc)

Now, for my next step, I am trying to add a code that gives users the option for showing the report on the screen or creating a text file.

What function do I need to do this?


texts: Batching field calculations and feature classes in feature dataset using ArcPy?

texts:
I'm working on writing a script that will calculate facility ID's for feature classes in a feature dataset. Is there a way to do this without making feature layers and doing select/field calculate operations on all of them?

Basically, I'm wanting to iterate through the feature dataset and create a list of feature classes that can then have the field calc ran on each. The problem is that each feature class has a different expression for field calculation, and I need to select only null values, and the select layer by attribute function requires a feature layer.


texts: How to convert polygon into hexagonal grid?

texts: I'm trying to convert a polygon shapefile of the UK into a hex grid. Using MMQGIS I've produced the grid, then I tried the intersect tool which gave me something close to what I want. However, the irregular edge of the UK borders has cut some of the hexagons. Is there a way of keeping all of the hexagons whole and producing an 'approximation' of the original polygon.


texts: QGIS project on QGIS server -- cannot load wms layer

texts:
I've created a qgis project on QGIS 1.8, copied it onto the QGIS Server, and am trying to load the project as a WMS Layer back into QGIS. There is one vector layer in the project. Paths have been set to absolute and the CRS for both the layer and project is EPSG 27700.

In the "Add Layer(s) from a Server" dialog box, the project is showing, but not the vector layer. Clicking on the project, at the bottom of the dialog box it says "No Common CRS for selected layers".

Any help is very welcome -- I have searched the forums and tried the solutions suggested QGIS Server Capabilities Response doesn't contain any Layers and QGIS Server Capabilities Response doesn't contain any Layers (Fedora 15 64-bit) but still cannot get it to work.

Many thanks!


texts: Restart IIS Admin Service or Recycle Default App Pool?

texts: A weekly maintenance process requires stopping ArcGIS Server and ArcSDE services. After restarting all the services, is it better to restart IIS Admin Service or recycle Default App Pool? I currently run a python script to do the later, with inconsistent results.  Sometimes it takes 2 or 4 attempts, or stopping and starting AGS and ArcSDE services again, to get the web services working.  I'd like a clear procedure that will work each time.


texts:
I have to generate LONGEST FLOW PATH(polyline) within watershed. What is the procedure to generate this ? After generating longest flow length , I will divide polygon with polyline.
I have :

watershed layer(polygon)
stream layer(polyline)
Digital Elevation Model (DEM) layer

I am using .NET and python in ArcGIS 10. I have tried "FlowLength()"  http://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/flow-length.htm
method , but it gives raster output. I tried to convert this raster dataset to polyline. this was my procedure.

I have attached an image to simplify this problem.


texts: Compare multiple rasters with consideration to order over time

texts:
I have multiple land cover rasters for each year (i.e. 2002, 2003, 2004...etc.) and each of them have five land cover classes. I would like to set up evaluation rules in ModelBuilder.

An example of a rule would be: can't have class 2 or 3 and then have class 1. To clarify, if 2002 had a value of 1, and 2003 had a value of 2, and then 2004 had 1 again, this would be flagged.

I tried the combine tool, to have all the years in one single raster, and then evaluate across the fields.

Is there a way to achieve the above scenario?

The important distinction here is that anytime a sequence like 1-2-1 or 1-3-1 occurs, no matter the years, I'd like to identify it.


texts: How to prepare so many raster data into one ascii format

texts: I have a number of spatial data on which I will like to perform an econometric analysis. The data include land use data, soil data, population data and so on which are all in raster format. I need to convert them into stata files to run the econometric model. I am confused on how to go about this. I need to prepare all the data in one ASCII text file for me to convert to stata files using another software. I need to prepare all the data as one sample, A guide on how to prepare the sample is in the image below. I do not know how to go about this.


texts: Given 2 LINESTRINGs with common point, how to extract common point?

texts:
I have 2 LINESTRING geometries. I know both share a common point. Using PostGIS for Postgres, what is the cleanest way to extract the common point?

Example with illustration:

How do I extract the intersection of the green and blue line, assuming the points at the crossing are equal?


texts: Is there a way to reduce tile processing from Openlayers 3 when user pan & zoom a lot?

texts:
I have 2 overlays with local tiles (xyz) on my map.
When the user zoom in / out and pan the map intentionally a lot, the app becomes very slow because OL tries to render all the tiles it was requested by the pan & zoom, but they are no longer needed on the view because the user zoomed out from there.

Is there a way to reduce this load? Will it be faster if instead of xyz layer I'll use geojson layer with styles?


texts:
I have a shp file that consists of the road system of Mexico City. I would like to count the number of road intersections (crossings) and blocks (which I would define as a segment starting from one intersection and ending in another). This within a buffer of a certain size: a circle centered in a polygon within another shp file.

I am mainly using ArcGIS but am quite familiar with R, Python and QGIS if the solution is simpler on any of those. I realize their are a few similar questions around here but havent been able to use any of the information provided there.


texts: How to find the service pack level in ArcGIS Addin?

texts: I can query the current version from RuntimeInfo.Version. But I would also like to know what service pack is also installed. Do I need to just check the registery key HKEY_LOCAL_MACHINE\SOFTWARE\ESRI\ArcGIS\RealVersion. Or is there a proper ArcObject that allows such a thing?


texts: Geoserver, PostGIS and Postgres setup better on which web server?

texts: Is it better to have Jetty or Apache as a Webserver for GIS stack?


texts: Performing Multicriteria Analysis using QGIS?

texts:
I have to do a multi-criteria analysis to answer the question : "which is the best lot to develop".

A few of the criterias are :

distance of the nearest bus stop (point layer with bus stops)
distance of the nearest shop (point layer with shops)
what is the flood danger (polygon layer, with danger grade attribute from 1 to 4)
is the lot in a nature-protection area (polygon layer)
is the owner already planning something on his lot (manual entered information in the lot's attributes)
and so on...

I thought I'd give it a try with QGIS, and here's how I've done :

add the following columns in my lots layer attributes table :

"analysis_BUS"
"analysis_SHOPS"
"analysis_FLOOD"
"analysis_PROJECT"
"..."
"analysis_MEAN"

Convert my lots layer to points using "polygons to centroids"
Run the "distance matrix" tool
Open the CSV to run an operation in excel (bus stop grade is 1.0 if nearer than 200m, and 0.0 if more than 750m, but I coulnt find the MIN() function in QGIS)
Join the resulting CSV back in QGIS
Repeat the same for shops
Run the "point in polygon" tool to select all the points in nature-protection area
Set 0.0 to all selected points
Repeat for other "in ... area" criterias
Run the "spatial join" tool to merge flood danger area information
Run a calculation using the column calculator to have the mean grade (using determined factors for each criteria)
Once all that done, re add the BUILDING LOTS shapefile once for each criteria
For each criteria, join the converted layer (that one with the centroids)  on the LOT id
Set the display to a gradient from red to green according to the corresponding criteria attribute and the mean grade attribute

Now, after a good 2 days of work, i now have all my criterias displaying in green if good choice for building, and red if bad choice, and I have my synthesis which aggregates all my criterias in one beautiful red-green map. (and I also have a huge mess in my "shapefiles" folder)

Now the problem.

What if :

i'd like to try the same analysis with another bus network scenario ?
i receive an updated lots shapefile (with, let's say, 13 modification in all the 13000 lots)
i'd like to test different weights for my criterias ?

Do I have to start all over again ?

Am I correctly using the wrong tool, or am I using the correct tool wrong ?

Would it be easier with a commercial GIS software ?

I see what answerers/commenters mean, and i didn't really think of using rasters.

However, the main question was more about the ability to try different scenarios or update the base data without having to restart all the process from scratch.

It seems that your suggestions are not much more flexible than what I suggested (even maybe more complex) since you have news steps :
- (for each criteria) rasterization.
- (in the end) sampling (quite complex if you want to include partial overlaps)

That Sextante Model builder seems awesome; in fact I was exactly thinking at something like that when posting my last comment.

I've used Grasshopper3D quite alot (it has nothing to do with GIS software) which is a great plugin for the Rhino3D modeler and which uses the same concept of node graph workflow construction. (example : http://designreform.net/2009/07/rhino-grasshopper-parametric-truss )

This seems so well adapted to a lot of GIS data analysis that I'd love to see a GIS software really built around such a node graph tool.

I'm looking forward to try Sextante Modeler and let you know how it worked out.
I wish i had found about it by myself by googling it, but i didn't know the keyword "model builder".


texts: Change Bootstrap Switch Depending on Leaflet Layer Zoom

texts:
I have this JSFiddle that as you can see has two tile layers. I have set the minZoom of the 'Example Layer' to 5 so that it only loads once you have zoomed in that far.

The issue is, using Bootstrap Switch I cannot get the switch for the 'Example Layer' to be disabled until the user zooms in to zoom 5 in this case and then the layer would be visible.

The ideal solution would add a tooltip to show on to the switch whilst it is disabled that advises the user to zoom in further to view the layer.


texts: Converting polyline to fixed width polygon using QGIS?

texts: I have a shapefile of polylines for trails, and I need to convert them all to be polygons that are 8 ft wide (4 ft from the center of each side of the polyline). I tried a following a few other posts that I thought were doing the same thing, but I always just ended up a big blob of color.


texts: How to set the coordinates for a polygon in openlayers?

texts:
I am very new to OpenLayers and JavaScript, and I have the following problem.

I have an .csv table representing the coordinates of points I want to visualize them on a map using OpenLayers.

I have found the following example on the OpenLayers page,

https://openlayers.org/en/latest/examples/polygon-styles.html

However, I couldn't understand the representation of the coordinates there.

More specifically, I didn't know what does the e mean in the coordinate
[-5e6, 6e6] for example.

However, I tried to plot a simple square on my map using my coordinates, but it just giving me a point, in the center of the map, as following:

https://jsfiddle.net/api/post/library/pure/#&togetherjs=bD5bfPm7vz

So I don't know what's exactly is wrong, and what should I change?

I think it's all in the way the coordinates are written, but not sure though.

And this is my code:

var styles = [
  new ol.style.Style({
    stroke: new ol.style.Stroke({
      color: 'blue',
      width: 3
    }),
    fill: new ol.style.Fill({
      color: 'rgba(0, 0, 255, 0.1)'
    })
  }),
  new ol.style.Style({
    image: new ol.style.Circle({
      radius: 5,
      fill: new ol.style.Fill({
        color: 'orange'
      })
    }),
    geometry: function(feature) {
      // return the coordinates of the first ring of the polygon
      var coordinates = feature.getGeometry().getCoordinates()[0];
      return new ol.geom.MultiPoint(coordinates);
    }
  })
];

var geojsonObject = {
  'type': 'FeatureCollection',
  'crs': {
    'type': 'name',
    'properties': {
      'name': 'EPSG:3857'
    }
  },
  'features': [{
    'type': 'Feature',
    'geometry': {
      'type': 'Polygon',
      'coordinates': [[[-5e6, 6e6], [-5e6, 8e6], [-3e6, 8e6],
          [-3e6, 6e6], [-5e6, 6e6]]]
    }
  }]
};

var source = new ol.source.Vector({
  features: (new ol.format.GeoJSON()).readFeatures(geojsonObject)
});

var layer = new ol.layer.Vector({
  source: source,
  style: styles
});

    var basic = new ol.layer.Tile({
      source: new ol.source.OSM()
    });

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


texts:
Moran’s I, a measure of spatial autocorrelation, is not a particularly robust statistic (it can be sensitive to skewed distributions of the spatial data attributes).

What are some more robust techniques for measuring spatial autocorrelation? I’m particularly interested in solutions that are readily available/implementable in a scripting language like R. If solutions apply to unique circumstances/data distributions, please specify those in your answer.

EDIT: I’m expanding the question with a few examples (in response to comments/answers to the original question)

It’s been suggested that permutation techniques (where a Moran’s I sampling distribution is generated using a Monte Carlo procedure) offers a robust solution. My understanding is that such test eliminates the need to make any assumptions about the Moran’s I distribution (given that the test statistic can be influenced by the spatial structure of the dataset) but, I fail to see how the permutation technique corrects for non-normally distributed attribute data. I offer two examples: one that demonstrates the influence of skewed data on local Moran’s I statistic, the other on global Moran’s I-–even under permutation tests.

I'll use Zhang et al.'s(2008) analyses as the first example. In their paper, they show attribute data distribution's influence on the local Moran’s I using permutation tests (9999 simulations).  I’ve reproduced the authors’ hotspot results for lead (Pb) concentrations (at 5% confidence level) using the original data (left panel) and a log transformation of that same data (right panel) in GeoDa. Boxplots of the original and log-transformed Pb concentrations are also presented. Here, the number of significant hot spots nearly doubles when the data are transformed; this example shows that the local statistic is sensitive to attribute data distribution--even when using Monte Carlo techniques!

The second example (simulated data) demonstrates the influence skewed data can have on the global Moran’s I, even when using permutation tests.  An example, in R, follows:

library(spdep)
library(maptools)
NC <- readShapePoly(system.file("etc/shapes/sids.shp", package="spdep")[1],ID="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
rn <- sapply(slot(NC, "polygons"), function(x) slot(x, "ID"))
NB <- read.gal(system.file("etc/weights/ncCR85.gal", package="spdep")[1], region.id=rn)
n  <- length(NB)
set.seed(4956)
x.norm <- rnorm(n)
rho    <- 0.3          # autoregressive parameter
W      <- nb2listw(NB) # Generate spatial weights
# Generate autocorrelated datasets (one normally distributed the other skewed)
x.norm.auto <- invIrW(W, rho) %*% x.norm # Generate autocorrelated values
x.skew.auto <- exp(x.norm.auto) # Transform orginal data to create a 'skewed' version
# Run permutation tests
MCI.norm <- moran.mc(x.norm.auto, listw=W, nsim=9999)
MCI.skew <- moran.mc(x.skew.auto, listw=W, nsim=9999)
# Display p-values
MCI.norm$p.value;MCI.skew$p.value

Note the difference in P-values. The skewed data indicates that there is no clustering at a 5% significance level  (p=0.167) whereas the normally distributed data indicates that there is (p=0.013).

Chaosheng Zhang, Lin Luo, Weilin Xu, Valerie Ledwith, Use of local Moran's I and GIS to identify pollution hotspots of Pb in urban soils of Galway, Ireland, Science of The Total Environment, Volume 398, Issues 1–3, 15 July 2008, Pages 212-221


texts: Tracing Boundary Around a Polygon Shapefile in QGIS

texts:
I'm looking to create a boundary around this shapefile of street blocks to have show the city limits in QGIS.

I tried this solution, (Settings > Snapping Options and enabling the "Avoid Int." checkbox, then adding a feature to the polygon via the editor tool) but couldn't get it to work. It is supposed to snap to the outermost parts of the polygon if I interpret it correctly...

Ideally, the solution would create a polygon that removes the spaces but keeps the outermost parts of the polygon. I believe there is a similar tool in ArcGIS called the Trace Tool if that helps. I'm hoping to get around tracing these by hand via the Snapping Tool.

Thanks everyone,

Zach


texts: How do I install the OpenGeo Suite on a web server?

texts:
I've been doing some research over the past few days to learn just how I can bring my parcel, right of way, easement etc... GIS data out to the field with me on my Android mobile devices with real-time GPS capabilities. I'm not in a position to spend hundreds of dollars a month for an ArcGIS subscription, so I'm trying to learn OpenGeo Suite instead. I've very limited experience with web servers, linux, etc... so I know I've got a lot of reading ahead of me, but I'm hoping someone here can put me on the right path or let me know if I'm just too far out of my league to do this without a lot of prior experience.

I do have a website and host running cpanel with a bunch of things already installed including Apache 2.2.26. If I'm understanding the OpenGeo instructions correctly, I need to install Apache Tomcat (or something similar) onto my web server in order to install OpenGeo's production WAR bundle? I see they also have installers for Ubuntu and Red Hat, however I have no idea if this is for a web server or desktop computer or either. Am I on the right track with this line of thinking or way off?


texts: Why do my Shapefiles not align correctly?

texts:
First, I'll summarize what I understand so far so please confirm or correct me.

When On-the-Fly is selected a layer is reprojected from that of the layer's CRS to that of the project CRS, if the project CRS has been defined    and the layer CRS has been recognized   (which apparently is not always the case);
The project CRS can be defined three ways: choose one from a menu, set it to a particular layer's CRS, let it default to a past CRS or a layer CRS.
The layer CRS may be recognized on import (or not).  If it's not (or is) the way to configure the the CRS is to save the layer as a shape file and select a CRS when prompted.  However, it sounds like if a .qrj file has already been generated this will fail somehow, so we must assure ourselves a .qrj file is not present.

How am I doing so far?

Second:  I'm using QGIS 1.8 on a Windows 7 OS.

Let's say two different layers have been saved each with CRS designated as EPSG:2272 (NAD83 / Pennsylvania South (ftUS)) and we start a new project.  We select On-the-Fly reprojection and set the project CRS to EPSG:2272.  We open both .shp files as layers.  They are both files of Philadelphia, (PhilaCopIncsJune302012, PhilaCensusBlocks2012) and the the .prj files respectively are:

 PROJCS["NAD83 / Pennsylvania South  (ftUS)",
 GEOGCS["NAD83",DATUM["North_American_Datum_1983",
 SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],
 TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],
 PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
 UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],
 AUTHORITY["EPSG","4269"]],
 PROJECTION["Lambert_Conformal_Conic_2SP"],
 PARAMETER["standard_parallel_1",40.96666666666667],
 PARAMETER["standard_parallel_2",39.93333333333333],
 PARAMETER["latitude_of_origin",39.33333333333334],
 PARAMETER["central_meridian",-77.75],
 PARAMETER["false_easting",1968500],
 PARAMETER["false_northing",0],UNIT["US survey
 foot",0.3048006096012192,AUTHORITY["EPSG","9003"]],
 AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","2272"]]

and

 PROJCS["NAD83 / Pennsylvania South
 (ftUS)",
 GEOGCS["NAD83",
 DATUM["North_American_Datum_1983",
 SPHEROID["GRS
 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],
 TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],
 PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
 UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],
 AUTHORITY["EPSG","4269"]],
 PROJECTION["Lambert_Conformal_Conic_2SP"],
 PARAMETER["standard_parallel_1",40.96666666666667],
 PARAMETER["standard_parallel_2",39.93333333333333],
 PARAMETER["latitude_of_origin",39.33333333333334],
 PARAMETER["central_meridian",-77.75],
 PARAMETER["false_easting",1968500],
 PARAMETER["false_northing",0],
 UNIT["US survey foot",0.3048006096012192,AUTHORITY["EPSG","9003"]],
 AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","2272"]]

The CRS found in the PhilaCopIncsJune302012 info document is WGS84 decimal degrees.

The CRS found in the PhilaCensusBlocks2012 XML document is

  NAD_1983_StatePlane_Pennsylvania_South_FIPS_3702_Feet

The two layers are far far away from each other.  Any ideas why?


texts: Error when export NDVI result to CSV table Google Earth Engine and "collection" is not defined in this scope

texts:
var Akdarya = ee.FeatureCollection("users/xodjayevshovkat/kasbi_nishan_wuas");

//Define a geometry.
Map.addLayer(Akdarya);

// Make the geometry a feature collection.

// Import S2 image collection.
var S2 = ee.ImageCollection('COPERNICUS/S2')
  // Filter collection.
  .filterDate('2018-05-01','2018-09-30')
  .filterBounds(Akdarya)
  .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 20);

// Define a cloud masking function.
var maskcloud1 = function(image) {
  var QA60 = image.select(['QA60']);
  return image.updateMask(QA60.lt(1))//.clip(aaaaaa);
};

// Create Clip collection function
var clipToShp = function(image){
  var clipped = image.clip(Akdarya)
  return clipped
}

// clip all images to this shp
var clippedCollection = S2.map(clipToShp)

//

// Define an index function (return only NDVI).
var NDVI = function(image) {
  return image.expression(
    '(NIR - RED) / (NIR + RED)',
    {
      'NIR': image.select('B8'),
      'RED': image.select('B4'),
    }).rename('NDVI').copyProperties(image, image.propertyNames());
};

// Apply the cloud mask function.
var NDVIcol = S2.map(maskcloud1)
 // Apply the NDVI transformation.
 .map(NDVI);

// Calculate the per-pixel median NDVI for all images in collection.
var NDVImed = NDVIcol.mean();

// Define vis params.
var NDVIvis = {min: -0.001, max: 0.9, palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
    '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'};

//

// create a collection of the NDVI
var ndviCollection = clippedCollection.map(NDVI)

// create a feature so we can iterate through it later
var fet = ee.Feature(Akdarya)

// add the NDVItimeseries to the SHP
var addPropNDVI = function(img, feature){

  var newf = ee.Feature(feature)
  var featureNDVI = ee.Number(img.reduceRegion(ee.Reducer.mean(),Akdarya))
  var theDate = img.get("DATE_ACQUIRED")//.format("Y-M-D");
  var ndviDate = ee.String("NDVI_").cat(theDate)
  return ee.Feature(newf.set(ndviDate,featureNDVI));
}

var featNDVI = ee.Feature(ndviCollection.iterate(addPropNDVI, fet));

//Make a collection of the features
var featureCollectionNDVI = ee.FeatureCollection([
  featNDVI])

// Export a .csv table of date, mean NDVI for Akdarya

Export.table.toDrive({
  collection: featureCollectionNDVI,
  description: "NDVI_Export",
  folder:"Time_Series",
  regions: Akdarya,
  reducer: ee.Reducer.mean(),
  seriesProperty: 'name_e',
  fileFormat:"CSV"
});


texts: Arcmap 10.3 crashing on launch using either existing MXD or just launching app

texts:
I recently upgraded from 10.0 to 10.3 and am having a LOT of crashes when trying to launch the app.  I have been using MXD's started in 10.0 but re-saved under 10.3 as well as ones started fresh in 10.3, both seem to experience crashes on launch.  I have also tried just launching the app with no document to start a blank MXD and also get crashes on launch.  Each time it is at the very end stages of the launching process, application windows have spawned, things seem fine, then right at the end it crashes.

I did find an old python 2.6 folder on my C: drive from 10.0 which I removed to see if that would help as I know 10.3 uses python 2.7, but no change.

It's getting very frustrating as launching the program has become tedious, with it crashing repeatedly, then mysteriously launching just fine, working for a bit, then I try and open another instance and I get crahses on the new instance, or I close the working instance down and try to restart the exact same MXD and I get non-stop crashing.

I am running Windows 7, have updated to the latest updates.


texts: Extracting Roads and Trails from SPOT Imagery

texts:
Is it possible to extract Roads and Trails vectors from SPOT Raster Tiff Image? Actually I have tried SPRING software to extract road vectors from this raster. Yes, It gives good result, But unfortunately it extracts road data only. It fails to extract trails. I have shared a sample data via Google Drive,FYI EPSG: 26913

Anyone had succeeded extracting low width trails/Roads? If yes, please guide me to achieve it.I want to extract the data as polyline or polygon data.(any Free software to do this task would be much helpful)

Link for Sample TIFF SPOT_SAMPLE TIFF


texts:
I am attempting my first use of the modeler from the Graphical Modeler from the Processing menu in QGIS. I made a simple model for cliping species data by sanctuary.It has two Vector Inputs, uses some different algorithms (buffer, Union, intersect etc).

Different questions arise:

I want the model to iterate over one layer. I would do this using the button "iterate over this layer" once you start the model. is there a better way to do so? for example it would be nice to iterate not over each Feature (geometry) of the layer but over each value of an Attribute (so if there's two Features with the same value in a given Attribute it would run it in one go and therefore only one Output file would be generated).

Second Question goes to file naming the Outputs:
I would like the Output files to be named after the above mentioned value of a chosen Attribute. is there a way to so?


texts:
I'd like to create some flow maps like the ones shown in Underdark's post:
https://anitagraser.com/2016/12/18/details-of-good-flow-maps/

I have a set of location points, say 50, and I'd like to create one "all points to one" map for each point, thus 50 maps in this example. And another set of "one point to all" map for each point, so 50 maps more.

Instead of creating all the extra line layers, I'd like to accomplish this using a combination of Atlas and Geometry generator symbols in QGIS.

However, I'm not finding many examples to see how or if I can do it.

So, is there a site where I can find an exhaustive set of examples of QGIS geometry generator?

This is a clip of the filtered output showing just the top 10 lines from each point. I created all the lines from a csv with the point coordinates using a WKT. Arrows and point symbols are created with the geometry generator. The rest is controlled by Atlas.


texts: How to select all parcels of a street? [ArcGIS Server]

texts: I have an Arcgis Server Map Service with a Street layer and a Parcel layer, I need to get the info of each parcel of a street that is selected by the user, does anyone knows how to do that? I'm a little bit lost. I'm using Arcgis Server 9.3, Visual Studio 2008 with Web ADF for .NET, thanks.


texts: How to make the thickness of traffic volume flow symbology proportional to traffic counts?

texts:
I can scale the bandwidth proportionally to the traffic volumes; however, I cannot show the directionality of the movements. Volumes for directional moves are stored in two different columns of a database. My workaround for this is to import the layer twice and use two different symbologies, one for each direction. Is it to possible to do this without having to import the same layer twice?

Note that I need to show the bandwidth starting from the link so I cannot use the same offset because it create blank space between the opposite directions (it seems to me that the offset start from the middle of the thckness).

Edit: in the uploaded images you can see that the offset depends on the width of the links, so I is not possible to find a single offset value matching for all the links


texts: How to combine rows in attribute table using ArcGIS for Desktop?

texts:
I have an attribute table with three columns "FID" "Shape" and "Class"

All records have the same shape (polygon). Each record has its own FID. I have three classes, "A", "B", or "C". I have twenty records for Class C and would like to combine them into one record.

I have tried doing a dissolve without specific anything in the Dissolve Fields area and the result was a shapefile with all the Classes combined into one, which I don't want. I also tried doing a dissolve while specifying Class in the Dissolve Fields area. This resulted in a shapefile with the same twenty records for Class C.

How do I get those twenty polygons into one record?


texts: Combining multiple overlapping rasters - retain maximum value?

texts:
I have 40 overlapping raster grids (derived from flood modelling software) that I want to combine in such a way so as to retain only the maximum value in any of the given rasters at each grid location, in a new raster grid.

Is there any application or command using the QGIS suite that would permit this without being unwieldy?


texts: How to store raster brick in memory, rather than on disk in R?

texts:
I am trying to run analysis on the yearly weather data files here: https://www.northwestknowledge.net/metdata/data/. I'm running into a problem where each yearly file for each variable requires 2.4GB memory if I read it in using ncdf4 package. I can use the raster packages brick function to create a pointer to the file on disk, but this is very slow.
How do I store a raster brick in my computer's memory, rather than just create a pointer? Also, does anyone have ideas on how to reduce the size of these files in R (they're about 100MB on disk)?


texts: Geoserver - How to set up a workspace when creating a new image mosaic using gsconfig?

texts:
I want to automate the process of creating workspaces and image mosaics by using gsconfig.

First I defined the list of new work spaces to be created.

f_workspace = list(set(lst_fold).difference(set(lst)))

Then, from the list (if not empty) I create workspaces.
Then I create new image mosaics. Everything works perfectly. However, I would like to assign another workspace to each new image mosaic? The name of the workspace and the folder name where the granules and the configuration files for IM are located are the same (ws). So, I want every image mosaic (ws) to be assigned a workspace (ws). Can somebody help me?

if not f_workspace:
    print("The list is empty!")
else:
    for ws in f_workspace:
        cat.create_workspace(ws,'http://localhost:8080/geoserver/'+ws)

        cat_wrksp = Catalog("http://localhost:8080/geoserver/rest", username="****", password="****")
        cat_wrksp.create_imagemosaic(ws, "/var/lib/tomcat8/webapps/geoserver/data/coverages/"+ws)


texts: I can not export the image to the drive in Google Earth Engine (GEE)

texts:
Whenever I try to export to the drive the image gives me the following error:

Error: Invalid argument: 'image' must be of type Image.

This is my code:

var mod11 = ee.ImageCollection("MODIS/MOD11A2");

// Specify name of country.
var country_names = ['Vietnam'];
// get fusion table with countries
var countries = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw');
// select the country from the fusion table
var vietnamCountry = countries.filter(ee.Filter.inList('Country', country_names)).geometry();

// set start date
var startdate = ee.Date.fromYMD(2000,1,1);
// set end date
var enddate = ee.Date.fromYMD(2017,12,31);

// filter the data collection
var temp =  mod11.filterDate(startdate, enddate)
  .sort('system:time_start', false)
  .filterBounds(vietnamCountry)
  .select("LST_Day_1km")

// convert LST to celcius
var toCelsius = function(image){
  var time = image.get('system:time_start')
  var celsius = image.multiply(0.02) // scale factor
                     .subtract(273.15) // from kelvin to C
                     .rename("Celcius")
                     .set('system:time_start',time)
  return celsius;
};

var CollectioninCelsius = temp.map(toCelsius)

var viz = {min:20.0, max:34, palette:"000b5e,96048f,4286f4,e2f442,960424"};
Map.addLayer(CollectioninCelsius.median().clip(table),viz,"temperature map")
Map.centerObject(vietnamCountry,5)

// Predefine the chart titles.
var title = {
  title: 'Temperature over time',
  hAxis: {title: 'Time'},
  vAxis: {title: 'Temperature (C)'},
};

// create the chart
var tempchart = ui.Chart.image.seriesByRegion(CollectioninCelsius,
                           vietnamCountry,
                           ee.Reducer.mean(),
                           'Celcius',
                            1000,
                           'system:time_start',
                           'PROJECT').setOptions(title);

// print the chart
print(tempchart);

Export.image.toDrive({
  image: mod11.select('Temp'),
  description: 'Temperatura',
  scale: 10,
  region: table,
  crs: 'EPSG:4326',
  folder: 'TCC',
  maxPixels: 10000000
});

https://code.earthengine.google.com/25c00e940fdcd572f8b0524a50b86d94

Table:
https://code.earthengine.google.com/?asset=users/eldileoliveira/Novapasta


texts:
I have published a layer on Geoserver by cascading a remote WMS. i have successfully got all the layers from remote wms and publish one of those layers. But I want to re-project the layer from its native srs wgs84 to declared srs EPSG 3011. But when i displayed this layer in my openlayers application, it did not displayed. Then i go back to geoserver and change the declared srs to wgs84 as well, this time the layer displayed successfully.
my question is, why the layer not re-projected ?. OR is there any other way to re-project the layer other than mentioned above?


texts: ArcGIS Point every 100 meters along poline

texts:
ArcObjects with ArcGis 10.1

I have constructed a polyline in WGS84 that has straight segments and curved segments. The line is heavily geodesically densified. Now I have to step through every 100 meters and perform additionally calculations at these points.

IConstructPoint.ConstructAlong looks promising but the distance is spat ref units which for WGS84 is degrees so that rules that out.

The best I have been able to do is as follows:
Step through every vertex in the line (which is way more than needed) and:

1. use IPolyline.QueryPointAndDistance to get the distance along the line in spat ref units.
2. Get a subcurve at that distance using IPolycurve2.GetSubcurve
3. Finally use IGeometryServer3.GetLengthsGeodesic of the subcurve to get the length in meters.

Seems like a lot of overhead. There must be a way to get a point every x meters along the line. Most of the functions I have looked at are based on length in spat ref units which is not what I am looking for.

Applicable code is as follows:

Dim pntAlongLine As IPoint = New Point
Dim distAlongLineToPnt As Double = 0, distFromPntToLine As Double = 0, IsToRightOfLine As Boolean
pPolyline.QueryPointAndDistance(esriSegmentExtension.esriNoExtension, pntOnLine, False, pntAlongLine, distAlongLineToPnt, distFromPntToLine, IsToRightOfLine)

'get a subcurve
Dim pPolyCurve As IPolycurve2 = CType(pPolyline, IPolycurve2)
Dim outCurve As ICurve = Nothing
pPolyCurve.GetSubcurve(0, distAlongLineToPnt, False, outCurve)

'create a polyline array with one element
Dim pPolylineArray As IPolylineArray = New PolylineArray
pPolylineArray.Add(pPolyline)

'create a linear meter unit
Dim pSpatRefFact As New SpatialReferenceEnvironment
Dim pLinearUnitMeters As ILinearUnit = CType(pSpatRefFact.CreateUnit(esriSRUnitType.esriSRUnit_Meter), ILinearUnit)

'setup a Geometry server
Dim pGeometryServer As IGeometryServer3 = CType(New GeometryServer, IGeometryServer3)
Dim pDoubleArray As IDoubleArray
pDoubleArray = pGeometryServer.GetLengthsGeodesic(pSpatRef, pPolylineArray, pLinearUnitMeters)

'the reference array will have geodesic distances for each polyline. in this case we passed in only one line.
distInMeters = pDoubleArray.Element(0)

This code is dependent on each vertex in the line but what I would really like is a point every 100 meters along the line.


texts: Using v.net.path to give route that spends less CO2?

texts:
i am trying to calculate the route that emits less CO2 (for a master thesis). I have a line layer with a field with CO2 and Lenght (and of course, with all the other fields like id speed etc):

The goal is to instead of calculating the less cost route, i want that the v.net.path gives me the route that spends less CO2, its more or less the same anyway.

What do i need to put in each field of v.net.path?

->"Input vector line layer (arc)" - just add the layer with the road network (and with the CO2 information per line).

->"Centers point layer (nodes)" - i created a nodes layer with v.net.nodes i assumed this is what the field requires.

->"Name of file containing start and end points" - ok, it has to be in a certain format, can you please write an example? I used this (json file): "id -932241 4859815 -931884 4858560" is it correct? (EPSG: 3857)

-> etc

ADVANCED PARAMETERS
-> "Arc type"
-> the next two choose "CO2"

etc

I think that the problem is in defining the start and end point.

I run the v.net.path it runs but gives me no results. Maybe i should use another coordinate system?

Is there any example out there (with everything WELL explained)?


texts: My ArcMap take about 3-4 minutes to start. It was good before with less than 10-20 sec start time. I have Dell Precision 6500 Laptop with 16 GB Ram with i7 processor. I could not figure out what is making it late? Can anybody help or had similar problem solved?


texts: Python GDAL: find out spikes and pits in DEM

texts:
I'm working with on a DEM (Geotiff) in Python GDAL, and I'd like to find out the spikes or pits in my terrain.
I have no idea how to proceed, does anyone have good algorithm to help me out?
Thanks,
eo


texts: I have n points, how do I find the closest one from a point p by road?

texts:
I'm currently working on a map with stores and clients. I want to find the closest store by road from  a given client, with the aim of simulating the creation of a new store.

Requirements:

precision is not an issue (low quality data is OK)
web based (API, some code server-side, or javascript)
free or cheap (no 100000$ GIS solutions...)

Of course I could use google maps api and try every store. But I have  a lot of stores and clients, and I would like to run the simulation several times, so it would sent far too many queries.


texts: Closing ArcMap getting started dialog using ArcObjects by C#?

texts:
I don't want the ArcMap "Getting Started" dialog to appear after arcmap starts up.

I have a login form in my Extension and after login I set a proper MXD for the user. So I don't want the "Getting Started" dialog to appear.

Has anyone done this task?


texts: 2018 USGS Geospatial PDF produces `ERROR 1: insufficient arguments for Marked Content`

texts:
Issue

When processing the 2018 Geospatial PDFs from the nationalMap.gov an error occurs that was not seen when processing the same quadrangle from 2015.

I should note:  there is no problem with gdal_translate, but an extra noisy ERROR that wasn't in previous versions of the GeospatialPDF

What — Processing USGS US Topos from 2018 with gdal_translate
Error — produces ERROR 1: insufficient arguments for Marked Content
Compare — 2015 USGS US Topo does not produce this same error
Compute time — the 2015 US Topo also takes less time complete in GDAL compared to 2018's version

Steps to reproduce

GDAL translate produces a TIFF as expected.

2015 GeospatialPDF

# wget
wget
https://prd-tnm.s3.amazonaws.com/StagedProducts/Maps/USTopo/PDF/CA/CA_Descanso_20150223_TM_geo.pdf

#gdal
gdal_translate CA_Descanso_20150223_TM_geo.pdf CA_Descanso_20150223_TM_geo.tif -co COMPRESS=LZW --config GDAL_PDF_DPI 200

2018 GeospatialPDF

# wget
wget
https://prd-tnm.s3.amazonaws.com/StagedProducts/Maps/USTopo/PDF/CA/CA_Descanso_20180904_TM_geo.pdf

#gdal
gdal_translate CA_Descanso_20180904_TM_geo.pdf CA_Descanso_20180904_TM_geo.tif -co COMPRESS=LZW --config GDAL_PDF_DPI 200

# Errors many times; similar to this
ERROR 1: Pos = 47918296, insufficient arguments for Marked Content
ERROR 1: Pos = 47918296, insufficient arguments for Marked Content
ERROR 1: Pos = 47918317, insufficient arguments for Marked Content
ERROR 1: Pos = 47918327, insufficient arguments for Marked Content
ERROR 1: Pos = 47918327, insufficient arguments for Marked Content
ERROR 1: Pos = 47918329, insufficient arguments for Marked Content
ERROR 1: Pos = 47918334, insufficient arguments for Marked Content
ERROR 1: Pos = 47918334, insufficient arguments for Marked Content
ERROR 1: Pos = 47918334, insufficient arguments for Marked Content
ERROR 1: Pos = 47918339, insufficient arguments for Marked Content
ERROR 1: Pos = 47918343, insufficient arguments for Marked Content
ERROR 1: Pos = 47920226, insufficient arguments for Marked Content
ERROR 1: Pos = 47920231, insufficient arguments for Marked Content
ERROR 1: Pos = 47920268, insufficient arguments for Marked Content
ERROR 1: Pos = 47920268, insufficient arguments for Marked Content
ERROR 1: Pos = 47920271, insufficient arguments for Marked Content
ERROR 1: Pos = 47946009, insufficient arguments for Marked Content
ERROR 1: Pos = 47946015, insufficient arguments for Marked Content
ERROR 1: Pos = 47947917, insufficient arguments for Marked Content
ERROR 1: Pos = 47954326, insufficient arguments for Marked Content


texts: How to save view so all maps are same when exporting maps?

texts:
I have been exporting a number of different projections for a study. However, the problem comes when I want to edit features or zoom in or out. When I then return to the Layout view, the map is not in the same position or orientation. This results in slightly different maps and when I put them together the difference in zoom and position of labels is quite obvious. I can of course attempt to match them by moving the map around and zooming in Layout view, but you don't get it exactly right (and the slight difference is so clear when flipping through maps in a presentation).

Is it possible to create a template, or something like it, in Layout view so that no matter what zooming or moving around you do in Data view, you can return to the exact same position to export an identical map?


texts: Union of polygons doesn't work 100 % correctly

texts:
I'm using the function gUnaryUnion to merge polygons in R. In one particuliar case this doesn't work correctly - some of the old polygon contours are still visible:

Can anybody help ?


texts: Using SQL Server Spatial without ArcSDE?

texts:
At work we are using ArcGIS 10.1 for fairly simple spatial data.  We don't use the more advanced items such as topology that are available in Arc.  We have a couple SQL Server databases for other kinds of data and have built some tools for taking that data and visualizing it in Arc.  There is some data that we use a File Geodatabase for that I want to put into SQL Server, so that it can be more closely tied to the other data we have, and so that when we have 2 or 3 people working in different parts of the database, it doesn't lock the whole database from editing (like the File Geodatabase we are using now does).

However, at the moment we don't have an Arc Server license, so although I can connect to the data, I can't edit it through Arc.  I have looked at other topics on here about SQL Server Spatial and Arc (this one for example), but I haven't found out how to successfully edit SQL Server data without a license.

My guess is the process would look something like this

create a local copy of the features I want to edit (maybe as shp files, or a file geodatabase)
edit as necessary
use external python libraries to delete the original features on the server
use shp2sql to upload the edited features

I feel like there must be tools that do this already, but maybe not for the price we want to pay.  An Arc Server license has been talked about here, but I have no idea when it will actually be implemented, so I don't want to hold my breath.  I know that there are other platforms that would be perfect for this, like QGIS, but that would be harder to convince people of than getting an Arc Server License.

Can someone point me in the right direction?


texts: Use rightclick to move/translate a feature in Openlayers 3

texts:
I'm rebuilding a working GoogleMaps Api v3 site to use Openlayers (3/4).

I want to make vector points (coming from a Postgis table) draggable by rightclicking them, so that they can be moved (translated) on the map to another position. The new position should then be updated in the database (by calling a php-script and passing the new position).

I want to use a rightclick instead of a singleclick or click to prevent users from changing positions accidentally while dragging around the map. As this could occur in the drag-features example. On our GoogleMaps site this works perfect and users feel fine with rightclicking to make features draggable.

I was looking around to find a solution, since Openlayers until now does not offer a rightclick event. The only good example for rightclicking that always appears is the ol3-contextmenu.

Is this really the only work around? I don't need and I don't want to use context menus.


texts: How to set up a Web Map Service with QGIS?

texts:
I am very new to QGIS environment.I wanted to work with creating Web Map Service using QGIS server.

Is it possible to create Web Map Service in QGIS server?


texts:
I am completely new to GIS and therefore have questions regarding how the GIS Web development should be made, because the entire concept still confuses me.

So, imagine I have a map. And by this, i mean the shape files for the layers and so on.

How i should i use this in my web application?

I have read about MapServer, OpenLayers, PostGIS, but can't quite understand how and why this all enters the question.

If i have my map, my files, shouldn't the problem be just on to represent it in the web page?


texts: ENVI shows colors of different RGB code the same in a bitmap image wrongly?

texts:
As you see in this picture which I've captured in MATLAB, there are two colors.

Red with the RGB code (255,0,0)
Brown with the RGB code (138,139,1)

I've tried to show the same picture in ENVI, but as you see the RGB codes (255,0,0) and (138,139,0) both appear as red in ENVI.

Does anyone know why this occurs and what should I do to fix it?


texts:
I am trying to wrap my ahead around the meaning of the M-value.
I read ESRI's dictionary definition here: https://support.esri.com/en/other-resources/gis-dictionary/term/m-value

and was not able to get a good grasp of the concept.

What does the M-Value in GIS data represent and can you provide an example of using it or how its useful?


texts: Is there an analog of PhotoShop's "Magic Wand Tool" in ArcGIS Desktop?

texts:
I want to interactively select rasters pixels based on it's values and modify them.

Here is example of what this "Magic Wand Tool" should do:


texts: Great circles path not curved (in R)? Conventions and convertions problems?

texts:
I am given the great circle equations of two paths. The first one is

  Longitude = 168 degree.

and the second one is:

  tan(Lattitude)=.11 cos(Longitude) + .64 sin(Longitude)

However, I am not sure I understand what these represent.

For example, when I then try to depict these two great circles on a map
(For the code below, I use as model the example found here):

library(CircStats)
library(geosphere)
library(maps)

great_circle_1 <- function( lattitude ){
    return( 168 )
}
great_circle_2 <- function( longitude ){
    return( atan( ( 0.11 * cos(longitude * pi / 180) + 0.64 * sin( longitude * pi / 180 )) ) * 180 / pi)
}
map("world", col="#f2f2f2", fill=TRUE, bg="white", lwd=0.05)
lat_c1_1 <- -50
lon_c1_1 <- great_circle_1( lat_c1_1 )
lat_c1_2 <- 90
lon_c1_2 <- great_circle_1( lat_c1_2 )
inter_c1 <- gcIntermediate(c(lon_c1_1, lat_c1_1), c(lon_c1_2, lat_c1_2), n=50, addStartEnd=TRUE)
lines(inter_c1)
lon_c2_1 <- 150
lat_c2_1 <- great_circle_2(lon_c2_1)
lon_c2_2 <- -50
lat_c2_2 <- great_circle_2(lon_c2_2)
inter_c2 <- gcIntermediate(c(lon_c2_1, lat_c2_1), c(lon_c2_2, lat_c2_2), n=50, addStartEnd=TRUE)
lines(inter_c2)

My problem is that contrary to the examples on the website (and what my intuition would dictate), the two grand circles appear as a flat lines (rather than curved ones). I wanted to make sure that this is not due to a comprehension error on my behalf over how to understand the great circle equation as formulated.

I am assuming spherical earth.

I gather that there are different conventions as to the coordinate systems. Since this is a question about conversions, for information, I use the longitude and latitudes as given on google map, for example the original Waterloo has coordinates (50.7167, 4.3833).


texts: How to Define projection and enter the Geographic Transformation parameters from SGL-GEO to WGS84 in QGIS?


texts:
I am slightly confused as to why when i am working in projects the CRS changes when I add new data, e.g. add Bing maps or try to georeference an image.

The CRS is set in the project and when I am prompted I always choose the same one, expect that it seems to end up on another one!

I know that using the same CRS is very important and I am therefore very confused as to why this keeps happening. How can I stop it from changing?


texts: Using GDAL, NetCDF4 to GeoTIFF wrongly rotates the output GeoTIFF file 90 degrees

texts:
I have written a code that transforms a NetCDF4 data (it can be downloaded from here) to GeoTIFF file using this link. This code works fine. But it seems that it rotates my data 90 or -90 degrees:
import sys
import numpy as np
from osgeo import gdal, osr, gdal_array
import xarray as xr

filepath = '/home/gpm/3B-DAY.MS.MRG.3IMERG.20170216-S000000-E235959.V06.nc4'
var_name = 'precipitationCal'
src_ds = gdal.Open(filepath)
if src_ds is None:
    print('Open failed')
    sys.exit()

# if len(src_ds.GetSubDatasets()) > 1:
# if there is more than one variable in the netCDF4
subdataset = 'NETCDF:"' + filepath + '":' + var_name
print('subdataset:', subdataset)
src_ds_sd = gdal.Open(subdataset)

print('opening dataset successful', src_ds_sd)

# getting info of the variable (subdataset)
ndv = src_ds_sd.GetRasterBand(1).GetNoDataValue()

xsize = src_ds_sd.RasterXSize
ysize = src_ds_sd.RasterYSize

geo_t = src_ds_sd.GetGeoTransform()
# geo_t = (geo_t[3], geo_t[5], geo_t[4], geo_t[0], geo_t[2], geo_t[1])
projection = osr.SpatialReference()
projection.ImportFromWkt(src_ds_sd.GetProjectionRef())

# close the subdataset and the whole dataset
# src_ds_sd = None
# src_ds = None

# read data using xrray
xr_ensemble = xr.open_dataset(filepath)
array = xr_ensemble[var_name]
array = np.ma.masked_array(array, mask=array==ndv, fill_value=ndv)
# array = np.transpose(array, (0, 2, 1))

final = (ndv, xsize, ysize, geo_t, projection, array)

datatype = gdal_array.NumericTypeCodeToGDALTypeCode(array.dtype)

if type(datatype) != np.int:
    if datatype.startswith('gdal.GDT_')==False:
        datatype = eval('gdal.GDT_' + datatype)

new_file_name = var_name +'33'+ '.tif'
zsize = array.shape[0]

# create a driver
driver = gdal.GetDriverByName('GTiff')
# set nans to the original no data value
array[np.isnan(array)] = ndv
# setup the dataset with zsize bands
dataset = driver.Create(new_file_name, xsize, ysize, zsize, datatype)
dataset.SetGeoTransform(geo_t)

spatial_reference = osr.SpatialReference()
spatial_reference.ImportFromEPSG(4326)

dataset.SetProjection(spatial_reference.ExportToWkt())
#%%
# write each slice of the array along the zsize
for i in range(zsize):
    dataset.GetRasterBand(i+1).WriteArray(array[i])
    dataset.GetRasterBand(i+1).SetNoDataValue(ndv)

dataset.FlushCache()

I am sure that I should change geo_t in some way. Currently the shape of the array is: (1, 3600, 1800)
3600 values for lon, 1800 values for lat
Here is the output geotiff file in QGIS:
. This is the output geotiff file's properties in QGIS:


texts: I have an openlayers map with more then 20000 features not in wms, only in vector layer. I have a control to select features, but sometimes it's really hard to select, specially the lines, and also the other feature because they are too close. So instead the regular select tool I want to be able to click on the map, and get all the close features then chose what to select. How I wanted to do it it's like this: when I click it's create a little circle and I check all the features that in the circle (even little part from the feature) I didn't really find a function to do it in the api, it's possible to do it? Or someone have other solution? Thanks


texts: Issue while performing merge/mosaic of multiple TIFF files using Rioxarray Python

texts:
I am trying the following code to mosaic/merge 10 image files (.tif files) using the following code:
import rioxarray
from rioxarray import merge
from rasterio.plot import show

items = [r'D:\LC08_L1TP_146046_20180425_20180502_01_T1.tif',
r'D:\LE07_L1TP_145045_20180426_20180523_01_T1.tif',
r'D:\LC08_L1TP_144046_20180427_20180502_01_T1.tif',
r'D:\LC08_L1TP_144045_20180427_20180502_01_T1.tif',
r'D:\LC08_L1TP_143046_20180420_20180502_01_T1.tif',
r'D:\LC08_L1TP_143047_20180420_20180502_01_T1.tif',
r'D:\LC08_L1TP_143045_20180420_20180502_01_T1.tif',
r'D:\LE07_L1TP_145046_20180426_20180523_01_T1.tif',
r'D:\LE07_L1TP_145045_20180426_20180523_01_T1.tif',
r'D:\LC08_L1TP_144047_20180427_20180502_01_T1.tif']

elements = []

for val in items:
    elements.append(rioxarray.open_rasterio(val))

merged = merge.merge_arrays(elements, nodata=0.0)

image = merged.values
show(image)

After running the above code I am getting the following plot:

The image shown above is not the correct mosaic/merge image as the left and right part of the image is swapped because when I am loading the input files in ArcMap I am getting the following:

How can I resolve this issue?


texts: qgis_version_2: add layer or group of layers always on the top of the legendInterface / Insert the layer at the top of the ToC

texts:
I have a code which is adding raster layeror group of raster layers, its working fine.

I want to know if there is a way to add layer or group of the layers always on the top of the  legend Interface?

root = self.iface.legendInterface()
if status: # True
    group = root.addGroup("Group")
    layer = self.iface.addRasterLayer(urlWithParams, "new_layer", "wms")
    root.moveLayer(layer,group)

else: # False
    self.iface.addRasterLayer(urlWithParams, "new_layer", "wms")


texts: Reclassifying raster/extracting pixel values conditionall over multipolygon shapefile

texts:
I have a shapefile containing 38 polygons i.e. 38 states of a country. This shapefile is overlaid on a raster. I need to extract/reclassify pixel above a certain value, specific to each polygon. For example, I need to extract the raster pixels> 120 for state/polygon 1, pixels> 189 for polygon 2 etc with the resulting raster being the extracted pixels with value 1 and everything else as NoData. Hence, it seems like I need to extract first and then reclassify.

I have the values, for extraction, saved as a data frame with a column containing names, matching the names of the states,which is stored as an attribute "Name" in the shapefile.

Any suggestion on how I could go about this? Should I extract the raster for each state into 38 separate rasters, then do reclassify() and then mosaic to make one raster i.e. the country?


texts: Is it possible to automatically set CRS of Rasters from PRJ file?

texts:
I am trying to learn my colleagues to display measured data (ASCII) in maps relatively easily. I can prepare project, color and symbology style for the point layer and background map.

I also found good way to get background map (PNG + world file) for offline use via MOBAC and OSM. MOBAC uses EPSG:3395 WGS 84 / World Mercator and I have the appropriate PRJ file. However, I found out, that QGIS does not use the PRJ for rasters and load them always using default CRS (in my case EPSG:4326).

Is it possible to make QGIS load the CRS settings from the PRJ (or other file) without need to set the CRS manually for every raster?

thanks


texts:
I don't know how to open up the log in Whitebox GAT so I can view and fix errors with running tools.

I tried to run the IDW tool for LiDAR data (inputting all necessary parameters), then an error message popped up saying "An error has occurred during operation. See log file for details." I I looked everywhere in the Whitebox options, under services, view, etc. but couldn't find a way to pull up the log


texts: Import multiple GML files into PostgreSQL 11

texts:
I have some 50 GML files in my downloads folder, each within a child folder such as:
Downloads\Adur\Land_Registry_Cadastral_Parcels.gml
Downloads\Abertawe_-_Swansea\Land_Registry_Cadastral_Parcels.gml
Downloads\Allendale\Land_Registry_Cadastral_Parcels.gml
etc...

I want to import all 50 GML files into a public schema as a batch command. I want to name each table as per the child folder (Adur, Abertawe_-_Swansea, Allendale etc) rather than the GML file itself. The DB has the PostGIS extension enabled. I'm using pgAdmin 4 as my visual SQL DB and expecting an idiot's guide on how to do this.


texts: Creating ArcMap Model with multiple iterators?

texts:
I'm using ArcMap version 10.4.1. I'm attempting to create a model which almost works perfectly but I can't figure out the last bit.

Ultimately I want to run the "extract by mask" tool multiple times. However, the input for the raster layer is in one folder and the input for the mask layer is in another folder. Each raster/mask combination represents one day; the data differs for each day. There are an equal number of raster and mask files in each folder.

This is the full model which I created. The model runs the whole way through without throwing any errors. When the model finishes I get different raster layers but all extracted to the same mask. Everything else works correctly- I get a folder of different masks, and I get another folder of different rasters. I want to extract each raster to the mask for that day.

The submodel on the top left creates the mask, while the submodel on the bottom left iterates through the folder which I'm populating the masks with. The submodel on the top right re-scales the raster layer. The submodel on the bottom right iterates through the folder where I put all the raster layers, and then runs the extract by mask tool. I put in pre-conditions (the dashed lines), so that the submodel on the bottom right, with the extract by mask tool, runs last.

The submodel on the top left, to create the cloud masks:

The submodel on the top right, to re-scale the rasters:

Submodel from the bottom left, to iterate through the folder of masks:

The submodel from the bottom right. This iterates through the folder of rasters and runs extract by mask. The white "input raster" bubble connects to the output from the model that iterates through the folder of masks.