texts: Setting up Generalized Linear Regression in ArcGIS Pro

texts:
I am analyzing a series of old canals and I am trying to make the case that these are tied to certain physical characteristics (soil, slope, etc) in the watershed. To check wether or not my hypothesis is correct I need to do some statistical tests. I used the Ripley's K-funtion tool in ArcGIS pro and did show indeed statistical significant clustering.
Now I am trying to use Generalized Linear Regression to prove the link with land. Since (to my understanding) the line-features are not supported in the GLR-tool and the canals starts in other landscape types than their ending ( wetter vs drier, more slope vs less slope, etc) I extracted the beginning and end of each point and added the different landscape variables to the point features. I generated several random points and merged it with the canal points. The random points got value 0 (no canal) and the canal points got value 1( yes canal) since some variance is necessary for the dependent variable. I ran binary GLR and Continuous (just out of curiosity) and the probability, VIF and Jarque-Bera stats are fine, however my adjusted R-squared is only 0,28.
I am assuming this is because with only the points instead of the lines, the link between landscape transitions and the canals is completely lost and the no-canal points are also present in the same landscape though not near the transitions like yes-canal points.
Any ideas on how to remedy this or evaluate my theory?


texts: export map state as JSON: Web_Map_as_JSON

texts:
Is there a way to export the current map state as a JSON object, mostly named Web_Map_as_JSON.
The JSON specification should be like this.

In other words, it would be convenient to have a function for the Map object, like myMap.Web_Map_as_JSON() that would provide the current state of the MAP.


texts: How to insert non-spatial data through WFS-T insert?

texts:
When i perform wfs-t insert operation using postgis database only spatial data and primary key gets inserted.

Can anyone point me to some examples in which attribute data are also inserted along with spatial data using openlayers?


texts: Calculating areas of polygons inside other polygons with GeoPandas

texts:
I have 2 geodataframes:

GDF1 is an area with overlapping polygons divided into 3 classes: (cleaning, infra and maintenance) (the black lines donates the overlapping polygons)
GDF2 consists of polygons describing the type of the area's
surface, like gras, sand, buildings, water etc (the colors donate the
type of area's surface)

I want to calculate: what kind of land_use_type (from GDF2) each polygon from GDF1 consists.

What i tried is :

GDF1 = gpd.GeoDataFrame.from_postgis('select class, geom_sm as geom
                                        From geo.GDF1',
                                        con=engine, geom_col='geom')

GDF2 = gpd.GeoDataFrame.from_postgis('SELECT land_use_type,wkb_geometry_sm as geom
                                      FROM geo.GDF2',
                                      con=engine, geom_col='geom')

I calculated the area of each land_use_type polygons:

GDF2["area"] = GDF2['geom'].area

And after that i merged on GDF1

merged = gpd.sjoin(GDF1, GDF2, how='left', op='intersects')

I got a new column with areas but I am not sure if that's the right calculation.


texts: Calculating/estimating carbon content in soil by using Remote Sensing

texts:
I am analysing sugar cane crops with QGIS and need to estimate the soil carbon content to evaluate the nutrition status of the crops.
The problem is that I don't have ground data so there is no possibility of doing any calibration. The only available data that I have is the shapefile with all the crops mapped as polygons and the NDVI and NDMI calculated previously.
I am wondering whether I can get this info through Remote Sensing (by using either Sentinel-1 SAR data or Sentinel-2).


texts: ArcPad not recording attribute info on Junos

texts: I have some Juno SBs and 3B which are running ArcPad version 10.2.  We run a custom Applet and these units will record points and let the user fill in the attribute information however the attribute data is not saved to the table, the fields are blank.  We are successfully running this Applet on dozens of other Junos with identical specs with no issues.  We have tried a hard reset and restoring factory settings, including uninstalling and reinstalling ArcPad, with no change. Any ideas on what would cause this issue that would not be fixed by a factory reset?


texts:
I am looking for a way to calculate the distance between different points in a linestring.
In practice I created a linestring from a X and Y coordinate by means of

st_linestring(asNumericMatrix(dataset[dataset$Individual==1, 5:6])

I achieved to get the angle between different points with

PolygonAngles(st_linestring(asNumericMatrix(dataset[dataset$Individual==1, 5:6]))

Now I am looking for similar way to calculate the distance between different points.

Any suggestion?


texts: grass r.watershed to larger rivers parallel lines to same stretch

texts:
I'm trying create stream channels with r.watershed grass64. However, it shows two or three parallel lines for the same stretch of river. The river is large in some parts, but is just one river. Does anyone know if some tools to fix this exist? I tried to use v.clean with some parameters, but was unsuccessful.

Another option could be using a polygon layer for the river for these large parts and after this convert to line. This polygon already exists and is available. Would be a sort of 'reconditioning', how it is in ArcGIS, but using polygon instead lines.


texts: ArcGIS Portal not working with LDAP users

texts:
I have set user identity store of ArcGIS Portal administration page like following.
{
  "type": "LDAP",
  "properties": {
    "userPassword": "secret",
    "isPasswordEncrypted": "false",
    "user": "uid=gisuser,ou=systemusers",
    "userFullnameAttribute": "displayName",
    "userGivenNameAttribute": "givenName",
    "userSurnameAttribute": "cn",
    "ldapURLForUsers": "ldaps://10.0.0.5:636/ou=users,ou=ags,dc=example,dc=com",
    "userEmailAttribute": "mail",
    "usernameAttribute": "uid",
    "caseSensitive": "false",
    "userSearchAttribute": "cn",
  }
}

After save the configurations, I tried to add members based on existing enterprise users in Portal. But no users found. The search result is empty whatever I searched.
Why does this happen?


texts: How to create a time series GIF web map?

texts: I found a web map on the New York Times showing early spring conditions throughout the US. I've seen a much more complicated time-slider widget that can be used with ArcGIS Online web maps, however, I really like the simple lines in this map. What technology was used to create this GIF web map?


texts: Creating summary table using ArcMap 10.2 Report Wizard?

texts:
I've been having a play and reading the help about a part of ArcMap I don't really use and that is the Reports. I thought I would teach myself a bit more about this by creating a simple report, for which I have spectacularly failed!

I have a point dataset, say with 100 points and there is a text field that gives each point a type. For arguments sake I have types A, B, C and D.

I've been trying to create a single page that would give a count for each type, something that the Summary Statistics tool can do.  The end result would be a nice looking report with a simple table like:

Type  Count
A      25
B      25
C      30
D      20

That is it! Nothing more complex. Have I misunderstood what a report is actually about? I've tried setting Group Levels, Summary Function to count and Summary Group properties and have not succeeded in creating a simple table.


texts:
I am looking to make a map from our online ArcGIS sites and host them on our own server. I got that all done, but I am trying to modify the code to my liking. There is a section in github (https://github.com/Esri/storymap-series) that goes over this which seems pretty self explanatory, sans one issue.

I can't seem to find where the map is declared. It says under application life cycle to look for the custom JavaScript code block at the far end of the index.html. The example they provide and the sites I got from our own ArcGIS does not have that code anywhere. Can't seem to find it in the .js files either.

Does anyone have any ideas?


texts: postGIS - Error: couldn't load library C:/Program Files/PostgreSQL/9.5/lib/postgis-2.2.dll: %1 is not a valid Win32 application

texts:
Environment - Windows10(x64), PostgreSQL 9.5(x64)

I installed GDAL, PROJ.4 and GEOS (via OSGeo4W )

Then I installed postGIS via Stack App Builder (I have tried both versions - x32 and x64). In both cases I have the same error - Error: couldn't load library C:/Program Files/PostgreSQL/9.5/lib/postgis-2.2.dll: %1 is not a valid Win32 application

I went to ../lib/ and postgis-2.2.dll is present there.

But when I am entering in psql shell - CREATE EXTENSION postgis; - I get this error.


texts:
I have a map with many feature classes on it. I want to use arcobjects to get to a selected layer (aka one of the feature classes) on the map, and change the fill color of all the polygons within that feature class. I don't want to add anymore features or make a new feature class, just change the colors of existing ones. I want to make this a user inputted color, so I want to use the color dialog in VS2013.

I've read that IGraphicsContainer is a useful interface to implement, but I don't know how to do that.

ILayerDescription layerdes = (ILayerDescription)featlayer;
layerdes.SelectionColor = colorDialog1.Color;

Implementing code with these lines gave me the error:
Cannot implicitly convert type 'System.Drawing.Color' to 'ESRI.ArcGIS.Display.IColor'


texts: ArcPy: Create NDVI with Raster Calculation and Map Algebra

texts:
My target is to create NDVI Rasters but I'm stuck with python loops, since I'm new to that, so sorry for a maybe basic question.

I have 2 folders, one with rasters (tif) with the red band, the other folder with the same rasters with the infrared band.

Folder A:
raster1.tif (red band)
raster2.tif (red band)
raster3.tif (red band)

Folder B:
raster1.tif (infrared band)
raster2.tif (infrared band)
raster3.tif (infrared band)

I want to create NDVI rasters with the following formula:
(NIR Band-Red Band)/(NIR Band+Red Band)

Of course, raster 1 from folder A has to match raster 1 from folder B.

I know I have to use map algebra (not Raster Calculator), but I'm already stuck with the loop, that the correct rasters are matching from folder A with folder B.

The following code doesn't produce the correct match. I tried it with a loop in a loop:

import arcpy
from arcpy.sa import*
arcpy.CheckOutExtension('Spatial')

arcpy.env.workspace = r"workspacetoredbandfolder"
liste_rotesband=arcpy.ListRasters()

arcpy.env.workspace = r"workspacetonirbandfolder"
liste_nirband=arcpy.ListRasters()

for rotesband in liste_rotesband:
    print rotesband
    for nirband in liste_nirband:
        print nirband
        break

but the result is:

raster 1- red band
raster 1- nir band
raster 2- red band
raster 1- nir band
raster 3- red band
raster 1- nir band

It should be:

raster 1- red band
raster 1- nir band
raster 2- red band
raster 2- nir band
raster 3- red band
raster 3- nir band

After this step, I was thinking about using the map Algebra to calculate my NDVI values.

Edit: I Tried it now with the following code:

import arcpy, os, math
arcpy.CheckOutExtension("Spatial")

workspaceRotesBand=r"workspacetoredbandfolder"
workspaceNIRBand=r"workspacetonirbandfolder"
OutputWS=r"workspacetooutput_ndvi_tiffs"

arcpy.env.workspace = workspaceRotesBand

rasterlist = arcpy.ListRasters()

for r in rasterlist:
    r1 = arcpy.sa.Raster(r)
    r2 = arcpy.sa.Raster(os.path.join(workspaceNIRBand,r))
    result = Float((Raster("r2")-Raster("r1"))/(Raster("r2")+Raster("r1")))
    result.save(os.path.join(OutputWS,r))

Now it seems that the linking works but I get the Error in the result line at the end: "'Float' is not defined".
If I remove the Float I get the Error "'Raster' is not defined".

Is my map algebra expression correct?


texts: OpenStreetMap Data Lookup Slow in PostGIS

texts:
I have imported OpenStreetMap data (Europe) into PostgreSQL 9.1 database with the help of osm2pgsql and filtered out everything expect administrative areas. I have GIST index for polygon (way) column. My planet_osm_polygon table has about 250 000 rows.

I want to make lookups to find all the administrative areas on which a set of coordinates is located. Currently I'm using SQL

SELECT osm_id,admin_level,boundary,name FROM planet_osm_polygon \
WHERE ST_Contains(ST_Transform(way,4326), \
ST_GeomFromText('POINT(val1 val2)', 4326)) ORDER BY admin_level;

But the query is too slow. How can I optimize the query and/or datastructure?

Some things that I have considered are:

Put bigger areas (countries) in an another table and connect them to the main table with foreign keys.
Change polygon datatype from geometry to geography (but then I have to convert all lookups because they are grades?).
Simplify the polygons with ST_SimplifyPreserveTopology().
Move index to a faster disk or memory if possible.


texts:
I am trying to connect with python code, QGIS and Spatialite, to do SQL query. I have read some tips and I tried with this code, but no add anything.

Any Idea?

uri=QgsDataSourceURI()
uri.setDatabase('C:\Users\pr\Desktop\GIS.sqlite')
schema=' '
table ='ambi'
geom_column= 'geom'
uri.setDataSource(schema, table, geom_column)
display_name=('prueba')
vlayer=QgsVectorLayer(uri.uri(), display_name, 'spatialite')

Finally I have put this:

uri = QgsDataSourceURI()
uri.setDatabase('C:/Desktop/GIS.sqlite')
schema = ''
table = 'PV'
geom_column = 'geom'
uri.setDataSource(schema, table, geom_column)
display_name = 'PV'
vlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
db = QSqlDatabase.addDatabase("QSQLITE");
db.setDatabaseName(uri.database())
db.open()
query=db.exec_("""select * from PV where pk=1 """)
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

But It adds the complete shape and I want only add, the PV with pk=1.
Could you tell me why it is?


texts:
Does an equivalent of ArcScripts exist?

I know of the the ArcGIS resource center and "Code Galleries" but does something that equates to ArcScripts where users share useful ArcGIS code exist? The new ArcGIS forums are not great for searching for code either.

It seems to me as through the ESRI replacement for ArcScripts are a poor replacement and I am wondering if something 3rd party exists.  I know we can still search ArcScripts but it does not appear to still be an active repository.


texts: FireCCI51 on Google Earth Engine: burnt area as shapefiles

texts:
How can I export my results as shapefiles (polygons) from the burnt area results (raster)?
I intend to export these shapefiles to ArcGIS for further anlysis on a specific national scale.
Link to example code:
https://code.earthengine.google.com/?scriptPath=Examples:Datasets/ESA_CCI_FireCCI_5_1


texts: Changing symbology of layer created using Select Tool using arcgis web app

texts:
We have an arcgis online app that allows users to select features with the option of 'create a layer'.  When the new layer is created from the selection it mimics the symbology of the original layer.

Is there a way for our customers to change the symbology of the layer they create?


texts:
This is my first time using DB Manager to store data.
The Import/Export them was easy.
My point is, when i try to run a QUERY it returns this error!

I´m using QGIS 2.18 , PostgreSQL 9.4


texts: How do the equations/algorithms behind buffering tools actually work?

texts: I've been wondering how the equations/algorithms behind buffering tools actually work?


texts:
I was trying to execute two tools (Mean Center Tool and Central Feature) when suddenly it gave me a error that I couldn't find the resolution too no matter what I tryed. The code is here (along with a link to the image of the error):

Traceback (most recent call last):
 File "c:\program files (x86)\arcgis\desktop10.2\ArcToolbox\Scripts\CentralFeature.py",       line 301, in <module>
 setupCentralFeature()
 File "c:\program files (x86)\arcgis\desktop10.2\ArcToolbox\Scripts\CentralFeature.py",   line 47, in setupCentralFeature
  ssdo = SSDO.SSDataObject(inputFC, templateFC = outputFC)
 File "c:\program files (x86)\arcgis\desktop10.2\ArcToolbox\Scripts\SSDataObject.py", line 288, in __init__
 self.distanceInfo = UTILS.DistanceInfo(self.spatialRef)
 File "c:\program files (x86)\arcgis\desktop10.2\ArcToolbox\Scripts\SSUtilities.py", line 2135, in __init__
 self.setInfo()
 File "c:\program files (x86)\arcgis\desktop10.2\ArcToolbox\Scripts\SSUtilities.py", line 2155, in setInfo
 info = distanceUnitInfo[self.name]
 KeyError: u'M'

 Failed to execute (CentralFeature).
 Failed at Thu Jul 31 16:10:13 2014 (Elapsed Time: 0.97 seconds)

Can someone explain what seems to be the problem here ?
If it's of good use, I am doing this on a point feature


texts: Intersecting raster with polygon grid in ArcGIS Desktop?

texts:
I'm trying to generate the fractional coverage of land-cover within polygon grid cells.
My land-cover is a 30m raster dataset (10008x8939 cells) and my grid is a fishnet shapefile (approximately 1km resolution; 120x111 cells).

What I normally do is create a fishnet corresponding to the raster, then sample the raster using the fishnet labels, and join the labels back to the grid polygon. This way I convert the raster into a polygon. Using the arctoolbox Raster to Polygon tool does not work so well for me because it does not give me a raster cell by cell representation of the data; even with "simplify polygons" unchecked.

With the raster as a polygon I can intersect it with the fishnet and do some calculations to get the fraction of land-cover types within each fishnet cell.

This worked really well until I got to this 30m dataset. Already ArcMap complains about a 2gb limit when trying to create the grid the raster will lay on.

Is there any other ideas?

A polygon grid that large will take ages to process.

The other answers related to this question seem to consider some sort of resampling and not relating the data to an arbitrary non-aligned grid.


texts: Qgis 2.2 error plugin manager/installer error

texts:
I get this error at startup

Ocurrió un error durante la ejecución del siguiente código:
import pyplugin_installer

Traceback (most recent call last):
  File "", line 1, in
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 217, in startPlugin
    msg = QCoreApplication.translate("Python", "%s due an error when calling its initGui() method" ) % errMsg
ValueError: unsupported format character ' ' (0x20) at index 2

Versión de Python:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

Versión de QGIS:
2.2.0-Valmiera 'Valmiera', c3a2817

Ruta de Python:
['C:/OSGEO4~1/apps/qgis/./python', u'C:/Users/Ing. Juan Bernnales/.qgis2/python', u'C:/Users/Ing. Juan Bernnales/.qgis2/python/plugins', 'C:/OSGEO4~1/apps/qgis/./python/plugins', 'C:\\OSGEO4~1\\bin\\python27.zip', 'C:\\OSGEO4~1\\apps\\Python27\\DLLs', 'C:\\OSGEO4~1\\apps\\Python27\\lib', 'C:\\OSGEO4~1\\apps\\Python27\\lib\\plat-win', 'C:\\OSGEO4~1\\apps\\Python27\\lib\\lib-tk', 'C:\\OSGEO4~1\\bin', 'C:\\OSGEO4~1\\apps\\Python27', 'C:\\OSGEO4~1\\apps\\Python27\\lib\\site-packages', 'C:\\Users\\Ing. Juan Bernnales\\.qgis2\\python\\plugins\\fTools\\tools']


texts: How can I make use of the QGIS authentication manager in PyQGIS processing scripts?

texts:
I have developped a set of QGIS processing scripts that rely on some background data stored in a postgreSQL database. During my development phase I had the database on my localost server and connected with hardcoded credentials in the QgsDataSourceUri connection. Now we have a common server and I want to enable my colleagues to use the scripts as well. I have migrated the database to the new server and each user who runs a script needs to authenticate with their own respective credentials. My first attempt to solve this is based on this great answer and works well so far:

uri = QgsDataSourceURI()
# assign this information before you query the QgsCredentials data store
uri.setConnection("HOST", "PORT", "DB_NAME", None, None)
connInfo = uri.connectionInfo()

(success, user, passwd) = QgsCredentials.instance().get(connInfo, None, None)

if success:
    uri.setPassword(passwd)
    uri.setUsername(user)
    uri.setDataSource("SCHEMA", "TABLE_NAME", "the_geom")
    layer = QgsVectorLayer(uri.uri(), "Layer", "postgres")

The drawback of this method is that the users need to enter their credentials each time they run a script, even when they have already authenticated themselves by entering their master password and already have access to the databases through the QGIS Explorer.

Is there a way to make use of the authentication manager in the processing scripts, so that the users only have to authenticate themselves once per session and not each time they run a script?

What security implications would this have?

I would be very greatful for any advice!


texts: Creating quarterly rasters from monthly rasters for multiple years with ArcPy?

texts:
I'm attempting to create quarterly rasters from monthly PRISM .bil files that I've already masked to my ROI. Essentially I'm trying to automate the process of using Raster Calculator to add 3 monthly rasters together (Jan + Feb + Mar = Q1) because I have 11 years of monthly rasters to combine into quarters.

I've already tried using an iterator in ModelBuilder, but I couldn't do what I need. I'm hoping this is possible using Python code though. My monthly rasters are organized into separate subfolders for each year with names like "07_0" for Jan 2007, "07_1" for Feb 2007, etc.

Is there python code that could iterate through each yearly subfolder and add each successive group of 3 rasters together into quarterly raters?

Here's my code so far:

import arcpy
from arcpy import env

arcpy.CheckOutExtension("Spatial")

# creat a variable for the folder with rasters
workspace = r'F:\Water Balance\Data\GIS_data\Rainfall\PRISM\Masked_rasters\2007'

# create a variable for output folder
ws_out = r'F:\Water Balance\Data\GIS_data\Rainfall\PRISM\Qtr_rasters'

# Set the environment workspace to the workspace folder
env.workspace = workspace
env.OverwriteOutput = True
# get a list of grids from the workspace
rasterList = arcpy.ListRasters()
# for loop to set groups of 3 rasters at quarters
for raster in rasterList:
    jan = rasterList[0]
    feb = rasterList[1]
    mar = rasterList[2]
    apr = rasterList[3]
    may = rasterList[4]
    jun = rasterList[5]
    jul = rasterList[6]
    aug = rasterList[7]
    sep = rasterList[8]
    oct = rasterList[9]
    nov = rasterList[10]
    dec = rasterList[11]

# use raster calculator sum the rasters in each quarter and make a new
raster
from arcpy.sa import *

# Use map algebra to add month raster into quarters
outRas = Raster(jan) + Raster(feb) + Raster(mar)
outRas.save("F:\Water Balance\Data\GIS_data\raster_test\2007_q1")


texts: Seeking inexpensive GIS software that can export TIN files to AutoCAD?

texts:
I'm an architect who would like to be able to download 3D terrain info and then bring it into AutoCAD 3D for use in renderings and animations using a couple of rendering/animation softwares called AccuRender and Lumion.

In the past, I was able to use an AutoCAD add-on software called EasySite (no longer available or supported) that could create TIN meshes from point files.

Is there an inexpensive way to download 3D terrain info that works in AutoCAD?


texts: I am using QGIS and HEC-RAS for my project. HEC-RAS needs a file for the river bathymetry. Since I do not have data for river bathymetry, I decided to use the Profile Tool in QGIS to create a river section. However, the Profile tool only allows to get a single cross section. What I want is the whole stretch of the river. Could anyone help me?


texts: I'm new to QGIS and stuck.  I have a shapefile of 1 ft contours for an area.  I want to load it into mapbox so I can turn it into a layer for Gaia GPS, but the shapefile is too big for mapbox and it gives an error.  I was told I could use QGIS to cut out a smaller piece from the contour map but I can't figure out how to do that.  I can load the contour map and see all the contours, but can't figure out how to select a portion of it and save to a new shapefile.


texts: Detailed changelog between minor QGIS version

texts:
How can I know in detail what fixes have been applied to a certain new QGIS version?

For example: what has been changed between QGIS 2.14.02 -> QGIS 2.14.03?

I know of the visual changelog, that gives a great overview of the large picture for 2.12 -> 2.14.

The answer is probably in the QGIS issues pages. But how to query it to know what fixes/improvements have been implemented between two certain versions?

@Joseph commented on the possibility to get the changelog from GitHub. Actually you get a good idea of the commits between versions by looking at the branch you are interested in, for example:
For QGIS 2.14 -> Then, viewing the raw file you can find "Release of 2.14.2", "Release of 2.14.1"... and in between the comments of the different commits that happened.


texts: SearchCursor() fails to read feature class?

texts:
I am very new to ArcPy and Python in general. I want write a script to read attribute table from a feature class and merge some data from a csv file into the table. Right now, I am just playing around with the table just to make sure everything is fine, so here is the code:

fc = 'ShuttleStop_nextbus'
cursor = arcpy.SearchCursor(fc, '*', """ 'NextBus_nu' = '27837'  """)
for row in cursor:
     print row

(I have set arcpy.env.workspace to current workspace)

So I inspected the content of the feature class with

len(arcpy.ListFields("ShuttleStop_nextbus"))

It return 13, which is exactly the number of fields. Then I read all the fields with

fields = arcpy.ListFields("ShuttleStop_nextbus")
for f in fields: print f

Then, it outputs something like this:

<geoprocessing describe field object object at 0x24546980>

<geoprocessing describe field object object at 0x24546DE8>

<geoprocessing describe field object object at 0x24546E78>

 .........
 .........

Turns out the fields are actually of type Describe.

Is there anyway that I can query the table if all the above occurs? I still want to use the "where" clause of SearchCursor(). Let me know if I need to clarify something.


texts: Converting osm file to shapefile (or data frame) in R

texts:
I was looking for a road shapefile for Australia today. I ended up going to Geofabrik. Open Street Map has data in .osm format. But they do not have one in .shp for Australia.

If possible, I want to convert the osm files to shapefiles using R. I have searched around, but I have not found solutions yet. Is there any way to to this? Alternatively, is it possible to read osm files into R and convert the data to data frame so that I can draw figures using ggplot2?

I have achieved the following graphic of NZ using two shapefiles with maptoolsand ggplot2. Ideally, I want to produce a similar map of Australia.


texts: Run function if item in QComboBox is selected

texts:
I got another issue writing a plugin to convert three shape files into an .obs input file. I created a UI (Qt) with a QComboBox to select a layer (that works fine). As soon as one of the layers in the QComboBox is selected, a second ComboBox should be filled with the fields of the layer to select from. To do that I'm trying to call a function, that is to be called, as soon, as a layer is selected in the first ComboBox. Unfortuntelly that does not work. Here is what I tried:

First I wrote a function to fill the second ComboBox with the field names:

  def build_Layer_Index_Changed():
    QMessageBox.information(None, "DEBUG", "Build Layer routine called. :)")
    canvas = self.iface.mapCanvas()
    dlg2 = ShowDialog2()
    build = canvas.layer(int(dlg2.Sel_Build_Layer.currentIndex()))
    if int(dlg2.Sel_Build_Layer.currentIndex()) >= 0:
      if not build.geomTypeDict == 2:   # Types: 0:"point",1:"line",2:"polygon", Buildings layer has to be of type 2
        QMessageBox.critical(self.iface.mainWindow(), ":-/ Selected buillding layer is not of type vector(polygon)! :-/" )
      else:
        build_fields = self.getFields(build)
        dlg2.Sel_buildhight_field.clear()
        dlg2.Sel_buildhight_field.addItems(build_fields)

"dlg2.Sel_buildhight_field" is the second QComboBox. I tried to connect it to the first ComboBox's (dlg2.Sel_Build_Layer) Signal by

QtCore.QObject.connect(dlg2.Sel_Build_Layer, QtCore.SIGNAL("currentIndexChanged"), self.build_Layer_Index_Changed())

but I get an error: "TypeError: build_Layer_Index_Changed() takes no arguments (1 given)". If I remove the brackets after self.build_Layer_Index_Changed, I get no eror, but the function is not called (as the MessageBox doesn't appear. I thought I did it exactly as given in the Example at http://techbase.kde.org/Development/Tutorials/Python_introduction_to_signals_and_slots, but obviously I'm doing something wrong…


texts: How to automate creation of Network Turn Features from line segments?

texts: I have a road network which I will add Turn features in. From the road network, I have created a model that cuts out the line segments where the turn features will be located. My problem is how to make a turn features of these line segments and get them to connect to the road network. I know how to draw in new Turn features, but there are too many turn features to draw here. There are approximately 28000 Turn features which should be included. Does anyone have a model, Python script or a solution on how to do this?


texts: Geowebcache bounding box to tile conversion is not working properly (EPSG:4326)

texts:
I am trying to make Geoserver and built in Geowebcache work together.
I have a page with google maps and I am rendering images from a published WMS layer using getTileUrl function.

Unfortunately I am getting wrong tiles when geowebcache is activated.

Here is a request which uses caching (caching is enabled for png8):

    http://myserver/geoserver/wms?REQUEST=GetMap&SERVICE=WMS&VERSION=1.1.1&TILED=TRUE&TRANSPARENT=TRUE&WIDTH=256&HEIGHT=256&FORMAT=image/png8&LAYERS=layername&SRS=EPSG:4326&bbox=28.7841796875,40.94671366508002,28.828125,40.979898069620155

response headers:

    Cache-Control:no-cache
    Content-Disposition:inline; filename=layername.png8
    Content-Type:image/png
    Date:Fri, 30 May 2014 08:19:17 GMT
    ETag:91591a0d685881053a9d4761f42086f8
    geowebcache-cache-result:HIT
    geowebcache-crs:EPSG:4326
    geowebcache-gridset:EPSG:4326
    geowebcache-layer:layername
    geowebcache-tile-bounds:28.7841796875,40.95703125,28.828125,41.0009765625
    geowebcache-tile-index:[4751, 2980, 12]
    Last-Modified:Thu, 29 May 2014 16:31:31 GMT
    Server:Apache-Coyote/1.1
    Transfer-Encoding:chunked

Here is a second request which hits directly to geoserver (caching is enabled for jpeg):

    http://myserver/geoserver/wms?REQUEST=GetMap&SERVICE=WMS&VERSION=1.1.1&TILED=TRUE&TRANSPARENT=TRUE&WIDTH=256&HEIGHT=256&FORMAT=image/jpeg&LAYERS=layername&SRS=EPSG:4326&bbox=28.7841796875,40.94671366508002,28.828125,40.979898069620155

response headers:

    Content-Disposition:inline; filename=layername.jpg
    Content-Type:image/jpeg
    Date:Fri, 30 May 2014 08:19:01 GMT
    geowebcache-cache-result:MISS
    geowebcache-miss-reason:no tile cache for requested format
    Server:Apache-Coyote/1.1
    Transfer-Encoding:chunked

Although I am requesting for the same bounding box, the first request which uses geowebcache returns wrong portion of the map while the second works fine.

Also, the first request finds a tile in the 12th zoom level while my map is actually in zoom level 13.

I guess, bbox -> tms tile mapping conversion is not working properly.
Anyone can propose suggestions on how I am doing wrong?

Note: my wms layer has native & declared srs:4326

Note 2: I am not using openlayers.

Edit:  I can check I am using the correct bounding box using the following url:

   http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/
   Google: (4751,3072)
   TMS: (4751,5119)
   Zoom:13

   == Tile Boundaries ==
   WGS84 datum (longitude/latitude):
   28.7841796875 40.94671366508002
   28.828125 40.979898069620155

But webgeocache finds [4751, 2980, 12].
In fact when I look at webgeocache directory tile x=4751 exists only in level 12 but not in level 13.

Edit 2: Example for topp:states layer
        Google: (301,382)
        TMS: (301,641)
        Zoom:10

    http://localhost:9090/geoserver/gwc/service/wms?REQUEST=GetMap&SERVICE=WMS&VERSION=1.1.1&TILED=TRUE&TRANSPARENT=TRUE&WIDTH=256&HEIGHT=256&FORMAT=image/png8&LAYERS=topp:states&SRS=EPSG:4326&bbox=-73.828125,41.50857729743935,-73.125,42.032974332441405

    response headers:
    Cache-Control:max-age=3600, must-revalidate
    Content-Disposition:inline; filename=geoserver-dispatch.image
    Content-Type:image/png
    Expires:Fri, 30 May 2014 10:47:17 GMT
    geowebcache-cache-result:HIT
    geowebcache-crs:EPSG:4326
    geowebcache-gridset:EPSG:4326
    geowebcache-tile-bounds:-73.828125,41.484375,-73.125,42.1875
    geowebcache-tile-index:[151, 187, 8]
    geowebcache-tile-index:[151, 187, 8]
    Last-Modified:Fri, 30 May 2014 09:43:38 GMT
    Server:Jetty(6.1.8)
    Transfer-Encoding:chunked

    Correct tile should be:
    http://localhost:9090/geoserver/wms?REQUEST=GetMap&SERVICE=WMS&VERSION=1.1.1&TILED=TRUE&TRANSPARENT=TRUE&WIDTH=256&HEIGHT=256&FORMAT=image/jpeg&LAYERS=topp:states&SRS=EPSG:4326&bbox=-73.828125,41.50857729743935,-73.125,42.032974332441405


texts: Unable to drag (re-order) layers within Table of Contents in ArcMap?

texts:
I added Raster data from a folder into ArcMap 10.1 and I cannot seem to drag any layer within the Table of Contents.

Can someone help?


texts: Using SQLite for multi-editors in ArcGIS?

texts:
I am wondering if it is possible to use a SQLite for multiple editors in ArcGIS Desktop?

ArcGIS 10.1
Small team, only five editors for now.
Shared network drive

Need a quick and solution due to governance and budget constraints while waiting for funding for a real SQL environment and IT team.


texts: Automating Atlas Generation With Different Coverage Layers

texts:
I'm currently using QGIS 3.10, and have around 60+ coverage layers that each should result in its own atlas pdf-file. All other settings in the print composer would be the same.

As I'm not too familiar with the QGIS API I'm having some trouble finding the right methods to use.

The easiest would be to write a short Python script to run on the interactive Python console in the QGIS GUI, and the logic I would like to implement would be something along the lines of this pseudocode:

make list of coverage layers
initiate print composer
load print composer template from .qpt-file (this is the same for all of the atlases)

for each LAYER in coverage layers
    set coverage layer to LAYER in print composer
    export atlas to pdf with filename LAYER.pdf

Any help is greatly appreciated


texts:
It contains two coordinate columns dist for distance in x direction and rel_alt for an altitude in y direction. The other two columns are the data points for each combination of those two coordinate values.

I am now trying to do a IDW interpolation on the ff column. This is what I tried:

#read the data
dat <- read.csv("data.csv")
#create the grid
x.range <- c(min(dat$dist), max(dat$dist))
y.range <- c(min(dat$rel_alt), max(dat$rel_alt))
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 1),
                   y = seq(from = y.range[1], to = y.range[2], by = 1))
#convert to SpatialPixels
coordinates(grd) <- ~x + y
gridded(grd) <- TRUE

neighbors = length(dat$ff)
beta = 2

idw_ff = gstat(formula = ff ~ 1,
                 data = dat,
                 nmax = neighbors,
                 set = list(idp = beta))

grid.ff <- predict(object = idw_ff, newdata = grd)

But I am getting an Error message saying:

Error in cbind(dataLst[[1]]$data@bbox[cbind(c(1, 1, 1, 1), c(1, 2, 1,  :
  trying to get slot "bbox" from an object (class "data.frame") that is not an S4 object

I don't get the Error, as the grd data is an S4 object. And I am not really sure if that what I am doing even would work properly?


texts:
I have been able to use PostGIS and Tiger successfully to geocode addresses that are more formally typed out, however when I try with something like

"Queens, NY" or "Brooklyn" it doesn't work.

SELECT * FROM geocode('brooklyn, ny');

I have seen this done, however, so I wanted to ask if anyone had some ideas on how I could do this? I would ideally like to get the first point or something of that sort to just get an approximate area of Brooklyn (maybe the middle, etc).


texts: Databases of bars, pubs and restaurants in United States and Canada?

texts:
I am mostly seeking Bars and Pubs but Restaurants are a very nice bonus. I need a database with zipcode, city, state, longitude, latitude, establishment name, establishment contact information, maybe a category with sub categories if it's restaurants as well. I found one with Oddity Software. But after closer inspection after buying it, the database was garbage. It had very long obtrusive category names and/or sub category names, that would make indexing it for any need difficult. Then it had so many entries that were not restaurants, pubs, or bars.. like auto repair centers, and a minimum of at least 80 categories/sub categories that made no sense to the type of database it claimed to be. All of which had multiple entries with that category/sub category that didn't apply.

So I see all these new sites popping up here and there trying to be like yelp or doing something completely different like I want to do. However I don't know where these people are getting their initial databases from. I thought Oddity Software would be my golden fix until I purchased it, found out what I found out about it, then they refused to discuss things further with me.

So I am wondering here today. Who knows of a database similar to that of which I seek?


texts: QGIS 2:14 - Border show-through on small Simple Markers

texts:
On small Simple Markers, 2.20000 mm or less, there is border show-through which I have been unable to eliminate with no pen or transparent border style.

This is shown in the group of large and small circles in the screen shot below.

An earlier post suggested using Font Marker or an SVG but these also have border show-through on small Markers

Any suggestions for solving this problem?


texts: Seeking purge function to remove layers from qgs file which are not present in Layers pane of QGIS?

texts:
I have some QGIS projects that I've been using for long time, adding and removing layers, updating reference to deleted or moved files etc.

Opening the qgs file with a text editor I see that it contains numerous references to different files, and this references are not showed on the Layer tree pane.

I have three different cases: references to existing file, references to missing files, references to database connections. I got aware of this when I opened my project on another machine, and QGIS reported a missing database connection to vector tables not present in the layer tree view.

Manually deleting references from the qgs file solved the problem, but I'm managing huge projects and it's near to impossible to manually check the project consistency. I've tried the "Remove empty layers" plugin, but the ones I'm trying to remove are not empty layers.

Is there any way to 'purge' qgis file?


texts:
I am trying to create an environmental sampling plan using an irregular shaped polygon. I need to divide an irregular shaped polygon into 30 equal area polygons. I have multiple sampling plans that I need to complete, so I have been searching for a way to automate this, rather than divide the polygon by hand.

Based on some research, I have found that creating a parcel fabric and using the parcel division tool will get me what I need. However, not only is my polygon irregularly shaped, but it also contains several donuts within the polygon. These areas represent islands and I don't want them included in the 30 equal area polygons.

I have tried to create this in two different ways and have yet to get it working. First, I created a parcel fabric excluding the donuts (there were no donut polygons - just holes of no data). When using the parcel division tool, I received the following error "A generated split line did not intersect with the source parcel" and the division was not completed. I assume this was because the split lines intersected the no data holes in the parcel.

I then created a parcel fabric including the donuts, each as individual parcels in the parcel fabric. I no longer am receiving the error, but the result I'm getting is just the large donut on its own and no division occurred.

If anyone has used the parcel division tool on polygons with donuts, how were you able to exclude the donuts from the parcel division?

To help visualize, below in purple are the donuts I want to exclude in the parcel division and the green is the area I want to divide.


texts: How to extract data from a polyline layer in a specified area(Arcgis)

texts:
I am trying to extract data from a polyline layer within set radii of preselected sites.  I have extracted other data from rasters using the Sample tool and rasters of the area of interest.  I would like to do something similar for polyline data.

Image showing clipped lines in green, original layer source in blue and 1km radius the layer was clipped by in black


texts:
I'd really appreciate some guidance with this. I need to create some map illustrated tiles to use in leaflet.js for a small web-app I'm building. So far though I've lost about five days to installing and trying to configure the GDAL package to work for me.

The things I have tried so far (all on mac OS 10.15.4):

download and installed GDAL2.4 package from http://www.kyngchaos.com/software/frameworks/
Once that was done, I scratched my head a lot, as I couldn't really find any documentation as to how exactly I'm supposed to use the package - sure it's installed, but what does that give me exactly?
After much searching and guessing, I found the script I wanted to run - gdal2tiles.py in my /Library/Frameworks/GDAL.Framework/Versions/2.4/Programs directory and from there was able to try running:
$ python gdal2tiles.py [args] but it returned an error:

  File "gdal2tiles.py", line 52, in
  from osgeo import gdal
  ModuleNotFoundError: No module named 'osgeo'

from there, I went down multiple rabbit holes

checking python versions (apparently I need to run 3.6 and I now am)
uninstalling/reinstalling all the things

But nothing fixes the problem.

Eventually I came across this post which suggests you can install GDAL2 via homebrew, which to me seems like it should be easy.. but no - whilst the process outlined was easy to follow, it left me with the same issue - I have GDAL installed, but no idea how to run the gdal2tiles.py script -I can still run the steps I outlined above and I get the same result. I'm guessing now I probably have multiple versions of GDAL installed which will bite me some other way.

Can anyone please help me with this? Essentially an idiots guide to installing and using GDAL would be perfect - all the info I can find during days of internet trawling seems to assume that I should know the answers already, but I'm not a Python dev, and have never run python packages before.

The main problems I would love solved would be:

The first issue - how do I install or make the osgeo module available to my gdal2tiles script?
If that's going to be too difficult, how can I access and use the gdal2tiles.py script through my homebrew installed version of GDAL2.

Some extra info

$ gdalinfo --version
returns GDAL 2.4.2, released 2019/06/28

If I open a python 3.6.5 shell, and try import gdal I get ModuleNotFoundError (someone suggested that this would test if gdal was "available" to my python install).

I've added a bunch of stuff to my .bash_profile as I've tried out different things and I have no idea what is helping or hindering anymore, so my .bash_profile now includes:

export PATH=/usr/local/bin:/usr/local/sbin:~/bin:$PATH
export PATH=/Library/Frameworks/GDAL.framework/Programs:$PATH
export PATH=/Library/Frameworks/GDAL.framework/Versions/2.4/Programs:$PATH
export PATH=/Users/smallmighty/Library/Python/3.7/bin:$PATH
export PATH=/Library/Frameworks/GDAL.framework/Programs:/Library/Framework/PROJ.framework/Programs:/Library/Frameworks/SQLite3.framework/Programs:/Library/Frameworks/UnixImageIO.framework/Programs:$PATH
export PYTHONPATH="/Library/Frameworks/GDAL.framework/Programs/"
export PYTHONPATH="/Library/Frameworks/GDAL.framework/Versions/2.4/Programs"


texts: Getting simple geoJson point layer into Leaflet map

texts:
I am having difficulty getting a geoJson layer with points to show up in my Leaflet mapp app. I want to add that the code is running on an online editor (Cloud9) and the geoJson file is local (on the same app as my app.js file). Here's what the folder looks like:

So far I have tried multiple ways to add the data to the map, using the Leaflet documentation. The first way was simply to use the following:

var myLayer = L.geoJSON().addTo(map);
myLayer.addData(./data/dtr_points.geojson);

I kept getting a parsing, unexpected token error in my eidtor. I then tried using the $getJson method, and even added styling to help it show up better, but nothing...

$.getJSON('./data/dtr_points.geojson', function (data) {
    // Assign the results to the geojsonData variable
    dtr_points = L.geoJson(data, {
      //inital css style of grid
      style: function(feature){
        return {
          color: "#000000",
          weight: 0.25,
          opacity: 0.5,
          fillOpacity: 10
        };
      }
    });
}).addTo(map);

I don't think there would be a cross origin issue with using an online editor. Can anyone suggest what I am doing wrong? Accoridning to Leaflet doc, points are a little different, but I am not sure why the above won't work.


texts: Convert lines to polygons, but closing the lines first (topographic lines not closed)

texts:
I have a topographic layer for a country. I want to select topographic lines of a certain altitude and transform the lines to polygons. All work fine except one thing: since the lines of the topographic map are not properly closed, the polygon are awful looking. See the image below.

As you can see, since the "topographic lines" are not joining the Country borders, it creates polygons that are closed from the perspective of the lines and not the country border. The blue line that I drew is the resulting polygon that I'd be looking for. So is there a way to tell QGIS to first close the topographic line to the border (small red line) and then to apply the line to polygon algorithm? Obviously, the small polygon would have the same treatment, but the image would be too cluttered if I drew another polygon to show what I'm looking for.


texts: QGIS 3.16: how to fix the position when open a form, Python Init function

texts:
When I open a form with "Identify Features", I would like to have it always at the same position (top left corner or in the center of the screen). Usually it moves every time I open it a little bit down and to the right and I have to get it back manualy to the start position.
I don't know Python, but I think I can do it somehow with the Python Init Code Configuration in the Layer properties. What code do I need to open the form on the top left corner or in the center of the screen?
I saw the question from user99 and question from kadeem, but I cannot bring this information to the Python Init Code Configuration from QGIS 3.16.
Which code do I need instead of the example code:
from qgis.PyQt.QtWidgets import QWidget
def my_form_open(dialog, layer, feature):
    geom = feature.geometry()
    control = dialog.findChild(QWidget, "MyLineEdit")


texts: Building raster indicating whether slope is positive or negative moving from start location?

texts:
I am using ArcGIS 10.1 (advanced user license), and I am wondering if there is any method to have ArcGIS calculate (for each cell of a slope raster in degrees) whether the slope is negative or positive, considering a movement from starting location outwards.

I found some threads (https://blogs.esri.com/esri/supportcenter/2012/03/09/calculating-the-slope-of-a-hiking-trail/; Calculating slope for points with downhill slopes rendered as negatives in ArcGIS) in which they were discussing a similar BUT NOT identical issue, and the solution was actually implying some coding I am not familiar with.

Is there any means to accomplish what I am after?


texts: Multiple Layers into Indiviual Shapefiles QGIS

texts:
I have extracted about 20 different tables out of an Oracle DB into a project now I want to save each one as a shapefile.
I know how to save an individual layer into a shapefile but is there a way to save them all at once.


texts: Are there any common rules for WMS clients to comsume the WMS in tiles?

texts:
Looking at OpenLayers and other clients. They requests small 256x256 tiles of the map. I am wondering if there is a common ruleset for these tiles to relax the server computation on raster layers.

Lets say a client need to show the map at BBOX, does it just generate a BBOX that fits exactly what it wants or does it adjust the box to fit some cache able bounding box that servers may implement?

My concern is that creating a WMS server won't scale well because clients can ask it to generate images as it fits. If the server site could just serve some cached tiles then it would scale much better.

I have seen there is something called WMTS, but our clients is asking for a WMS, so we can't just give them a WMTS.


texts: Load third-party APIs only when the corresponding map layer is requested?

texts:
OpenLayers allows multiple map layers to be hosted on a common platform. If the map application supports a wide variety of maps, say from Google, Yahoo and Microsoft, I find that upon initialization, the APIs from these map providers are loaded as well, regardless of whether I use them or not.

I am not sure if the calls to these third party javascripts count towards the usage quota of these map service providers (seems to be). Besides, there is also an additional delay in loading of the web page.

Having said that, I am looking for a way to put a placeholder in the layer switcher such that I can avoid declaring the third party map layer class and hence avoiding a call to the third party API script.


texts: Unexpected end of file when adding WFS on QGIS

texts:
I'm getting messages "Unexpected end of file" when adding a WFS on QGIS.
Should it be anything wrong with the URL or anything I could change?

Here is the URL which I tried to add and got the message.

http://sec.geodatenportal.sachsen-anhalt.de/wss/service/LVermGeo_INSPIRE_Raumbeobachtung_WFS/httpauth/licid-UUID_20140326-132415-051-81711-7868?Service=WFS&Request=GetCapabilities


texts: What do I need to create web based GIS application to show GTFS data

texts: I need to create several GIS based application and embed them on a web page for professional development. I need a cost effective solution and don't have big budget. Any help in right direction and to get started is appreciated (case study etc.).


texts:
Since the ability to store mxds inside a file gdb is still just an idea, I've been mulling over my options for using a "template" mxd within a model. My model uses 4 different mxds, but currently it accepts one at a time (i.e. I run the model 4 times for each feature class, using the different mxds). The model does various things such as saving a copy of the "template" as a new mxd, replacing data sources etc.

Up to now I have "Input Template mxd" as a model parameter. This is fine, if I'm the one running the model, so I know which mxd is the correct one. However, if someone else were to run the model, they could input any mxd they wished, which would inevitably lead to the model failing. The options I thought of were:

1. Trust that the user will enter the correct template mxd containing
   the correct layers: if only.
2. Force the user to enter a template mxd which has the correct name:
   won't work in this multiple template mxds scenario, as I would then have to
   duplicate the model 4 times to check for the 4 different names. Alternatively,
   if I wrote an if statement inside the model to check which template I'm receiving,
   it still would not potentially stop the user from renaming Template1.mxd to
   Template2.mxd, causing the wrong processes to run.
3. Store the templates inside the same file gdb with the model: as stated,
   cannot be done.

My current scenario has my model and template feature classes inside a gdb, inside a folder along with the four template mxds. The model automatically accesses each template mxd as needed, so the branch for each template runs consecutively everytime the model is executed. However, the user could still rename the mxds, or simply remove them from the folder...

My question is, is there a way to somehow store the template mxds within the model itself or in some other way so that the user cannot change it?


texts: Extracting cell size of Raster as value for buffering using ModelBuilder?

texts:
I am creating a model in ArcGIS ModelBuilder to buffer a line. I would like the buffer distance to be based on the cell size of a raster. The cell size may change depending upon which raster the user decides to use.

Is there a way to extract the value of cell size so that it can be used to populate the buffer distance field?


texts: Approach to controlling which layers display in a QgsComposerLegend

texts:
I've tried 2 approaches to control what shows up in print legends. Neither of them work which is troubling because other than pruning the layerRootTree in advance of composer.addItem(legend), there seems to be no way to display a layer in print without also displaying its symbol and name (as an entry) in the print legend's output.

    # APPROACH 1 - QgsComposerLegend(c).model()
    composerLegend = QgsComposerLegend(c)
    names = []
    excludes = ['100m', '10m', 'w001001']
    layer_set = []
    for node in self.getTreeLayers():
        layer = node.layer()
        if not layer.name() in excludes:
            layer_set.append(layer.id())
            names.append(layer.name())
    r.setLayerSet(layer_set)
    iface.mapCanvas().refresh()
    composerLegend.model().setLayerSet(layer_set)

    # APPROACH 2 - QgsComposerLegend(c).modelV2()
    composerLegend = QgsComposerLegend(c)
    composerLegend.modelV2().setRootGroup(QgsProject.instance().layerTreeRoot().findGroup('Land Uses'))
    composerLegend.synchronizeWithModel()
    composerLegend.updateLegend()
    c.addItem(composerLegend)

In the second approach the subgroups don't render properly in the legend and in the first all the layers show up regardless of the constraint assumably provided by setLayerSet.


texts:
Do you know of any good resources on geocoding algorithms?

I am particularly interested in address query parsing, including matching and weighting parts of the query, dealing with misspellings and variations, as well as in details about the physical data storage (e.g. schemas for direct relational database queries, approaches to data indexing etc.).

I have studied some documents about ArcGIS 10 geocoding, but they touch the actual implementation details just a little. Detailed documentation of other high-quality production implementations could also be helpful. The more technical the better. Theoretical algorithm papers are also great.

Thanks.


texts:
As you see on the attached image, I labelled contour lines on the map (in meters). Which function will allow the "halo" to take the background color of the map?

On the image, the halo is greenish which diminish the look of the map.

I use ArcMAP 10.


texts: Same colour palettes for two different maps in ggplot2

texts:
I am constructing two ion concentration maps in the UK using ggplot2. The dataset used for mapping are lzn.kriged1 and lzn.kriged2 which basically has the same structure:

> str(lzn.kriged1)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 9962 obs. of  2 variables:
  .. ..$ var1.pred: num [1:9962] 9.63 10.47 10.39 10.29 10.56 ...
  .. ..$ var1.var : num [1:9962] 10.59 2.16 2.1 2.11 4.79 ...
  ..@ coords.nrs : int [1:2] 1 2
  ..@ coords     : num [1:9962, 1:2] -6.35 -5.25 -5.19 -5.13 -5.72 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:9962] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "x1" "x2"
  ..@ bbox       : num [1:2, 1:2] -8.1 49.95 1.73 60.82
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x1" "x2"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr NA

I change the fill colour based on var1.pred with the following condition:

below 10 = color1
>10-20   = color2
>20-30   = color3
above 30 = color4

I use cut function to apply that condition:

lzn.kriged1 %>% as.data.frame %>% #lzn.kriged2 for MAP2
  ggplot(aes(x=x1, y=x2)) +
  geom_tile(aes(fill=cut(var1.pred,breaks=c(0,10,20,30,Inf),
            labels = c("Below 10", ">10-20", ">20-30", "Above 30")))) +
  coord_map() + theme(legend.title=element_blank())

which resulted in these maps:

In MAP 2, the maximum value of var1.pred is 19. That is why the legend only shows two levels.

My questions are:

1) How can I control the fill color for each level? For example:

below 10 = blue
>10-20   = green
>20-30   = yelow
above 30 = red

2) How can I make those two maps to have the same colors and the same legend?


texts: ScheduleViewer Google Maps JavaScript API error: ExpiredKeyMapError

texts:
I'm trying to validate a GTFS file with ScheduleViewer, but it gives me an error:

js?key=AIzaSyBLmKxIq4izZzlXKPeVoAiURcT6z6ynobg&sensor=false:51 Google Maps JavaScript API error: ExpiredKeyMapError
https://developers.google.com/maps/documentation/javascript/error-messages#expired-key-map-error
_.Oc @ js?key=AIzaSyBLmKxIq4izZzlXKPeVoAiURcT6z6ynobg&sensor=false:51

I changed the API key, but the error is not solved.


texts: Converting UTM to Lat/Long with a geoJSON file in pyproj

texts:
I am trying to convert a geoJSON file in UTM (EPSG 29902) to WGS84. I've seen other examples online which show how to use the transform function in pyproj; however all the ones I have seen are for a single set of coordinates which are typed into the console. For example:

wgs84 = Proj(init = 'epsg:4326')
IrishGrid = Proj(init = 'epsg:29902')
x1,y1 = 318909.45839999989, 251792.64200000092
x2,y2 = transform(IrishGrid,wgs84,x1,y1)
print (x2,y2)

Could someone please tell me how I would convert all of the coordinates in my GeoJSON file from UTM to lat/long? I have been able to read in the json file but am unsure of how to transform every coordinate. I would also like to be able to create a new GeoJSON file with the projected coordinates.

Here is my code so far:

from pyproj import Proj, transform

import json

with open('C:\Program Files (x86)\Anaconda3 (64-bit)\some_file.json') as f:
    data = json.load(f)

for feature in data['features']:

presume I use the transform function within my for loop but am unsure of the proper form.

The format of my json file is

"type": "FeatureCollection", "features": [ { "type": "Feature", "properties": { "REGION": "Southern Region", "REG_CODE": "03", "DIVISION": "Cork West", "DIV_CODE": "0319", "DISTRICT": "Bandon", "DIST_CODE": "4300A", "SUB_DIST": "Kinsale", "SUB_IRISH": "Cionn tSáile", "SUB_CODE": "4305B", "COUNTY_1": "Cork", "COUNTY_2": null, "GEOGID": "M4305B", "Male2011": 5765, "Female2011": 5963, "Total2011": 11728, "PPOcc2011": 4054, "Unocc2011": 1177, "Vacant2011": 1013, "HS2011": 5231, "PCVac2011": 19.4, "CREATEDBY": "Paul Creaner" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 154039.34449999966, 50090.991299999878 ], [ 154039.8311, 50105.332900000736 ], [ 154041.9757000003, 50130.335699999705 ],

followed by a list of coordinates to form a polygon for about 500 regions.


texts: ArcGIS Raster Calculator VS QGIS Raster Calculator

texts: Are there any major differences between ArcGIS Raster Calculator and QGIS Raster Calculator?


texts: Azimuth angle results are different for 'Field Calculator' and 'Measure Angle' on QGIS

texts:
I measured the angle between two lines with Measure Angle. The result is 39.531°. In attribute table I used this script for calculating azimuth with Field Calculator:

degrees(azimuth(start_point($geometry), end_point($geometry)))

The result of Field Calculator is 46.62005°. Why are there different results?

Coordinates:
Latitude1: 39.0390110  Longitude1: 33.8493433
Latitude2: 39.0413637 Longitude2: 33.8518330


texts: Attribute table missing after Refactoring fields in a layer

texts:
I wanted to add a new field to an attribute table for a county map I downloaded.  When I added it, I needed to change the location of the field to move it over next to the county name so I could keep track of it.

I downloaded table manager, but it came with a warning to use the refactor fields algorithm instead.

I used that, it worked. I saved the new table, added the new information. Everything was displaying correctly, I saved it, closed the program, and when I opened it back up, the refactored layer still appears, but there is no attribute table for it. It even saved all of the old categories I had specified.

Did I do something wrong?


texts: Calculating Polygon Width in ArcGIS Desktop?

texts:
Is there an easy to calculate the width of a polygon using ArcGIS 9.3.1 or 10?

I found some scripts for 9.1 but they will not work.  I can use width at widest point, narrowest point, or an average.  It does not need to be exact but needs to be somewhat automated because I have a couple hundred to calculate.


texts: How to render from Mapnik at a specific scale?

texts:
I'm trying to render a map to paper using Mapnik at a specific scale (such as 1:50k), but can't yet figure out how. So far I'm just using map.zoom_to_box:

map.zoom_to_box (mapnik.Envelope (101.5,3.2,102,3.3))

and play with MaxScaleDenominator and MinScaleDenominator in the styling rules..

But when done this way, I can only guess-work the map scale by adjusting the bounding box over and over, which is hard to obtain the exact 1:50k scale that I wanted.. I've also found this: map.scale_denominator() and map.scale(), but I think these are used only to display the map scale - not setting it..

Is there anyone here who can help me?


texts: How do I display geojson polygons on a map?

texts:
Here's a very basic question for you all.

I've got 17,000 polygons in a GeoDjango database. How do I create a slippy web map that can display these polygons?

In other words, do I:

load a GMap, get the boundaries, each time the map moves, make an Ajax call to the Django db, load the polygons...
make my own map tiles in advance (if so, how?)
something else?

I'd like users to be able to view the polygons on the map. I'd also like users to be able to search for a point, using GMaps geocoder or something similar, and load up the polygons within the borders of the map.

What client-side technology should I be using?

Thanks.


texts: Export.image.toDrive keeps failing for certain shapefiles in GEE

texts:
I am new to GIS. I had been trying to obtain village level images from kgis. But the GEE code to export images keeps failing with an error message of 'Internal error'. The code(Jupyter notebook) has worked fine for some district level shapefiles obtained from other website. Part of this code used for exporting:
# creating a task to export clipped image to drive
task = ee.batch.Export.image.toDrive(clip,
                                    description = 'image',
                                    folder = 'Landsat',
                                    region = geometry.geometry(),
                                    scale = 30,
                                    crs = 'EPSG:3857'
                                    )
task.start()

So, I am not really able to understand how to resolve this issue. Also the number of images obtained even after using filters are really high, which seems weird.


texts: Access to ArcObjects DLLs on a build server using Nuget

texts:
We're in the process of standing up a new build server--Jenkins on Windows Server 2012 if that matters--and are attempting to configure it to build a solution that uses ArcObjects (currently at 10.2.2, soon to upgrade to 10.4.1). In order to make the ArcObjects DLLs available to the build server we must install the ArcObjects SDK. In order to install the SDK we must install ArcGIS for Desktop, Server or Engine. Of course that software will never be used on this server so I'm looking for a way to swat our DLL fly without an ArcSledgeHammer.

I've considered circumventing the SDK on the server by tossing all ArcObjects DLLs into a nuget package on our internal repository. One down side of that would be duplicative references on development machines (which should still use the SDK). How great of a concern is that, and are there any other pitfalls to this solution?

We have at least five solutions that use ArcObjects:

an extension to ArcMap that does all sorts of editing & querying,
a server object extension,
a small application that is used to find & open ArcMap from other custom in-house applications,
another extension to ArcMap used for managing a specific dataset,
a standalone application that has the ability to export a shapefile (no spatial awareness in that application that I'm aware of).

The first one contains over 20 Esri.ArcGIS.Whatever references, and all the others likely reference subsets of that.  The first one and others (but not all of them) are currently being built by Hudson running on an old desktop machine that has ArcMap installed on it (because it's used now & then).


texts: Using Mapscript to create a mapfile from scratch

texts:
Following a few online tutorials I have seen, I am trying to create a mapfile from scratch using mapscript.  However, when I try to create the layer object, I get an error in the init of mapscript.py; NameError: global name 'args' is not defined.

Does anyone see something obvious?  When I use a mapfile that already exists on disk, everything works as expected.

        # create a new mapfile from scratch
        self.map = mapscript.mapObj()
        self.map.save('temp.map')
        self.map.setSize(256,256)
        self.map.setExtent(-180.0,-90.0,180.0,90.0)
        self.map.units= mapscript.MS_DD
        self.map.imagecolor.setRGB(255,255,255)

        # create a layer for the raster
        self.layer = mapscript.layerObj(self.map) ----> ERROR HAPPENS HERE
        self.layer.name= 'tilesource'
        self.layer.type= mapscript.MS_RASTER
        self.layer.status = mapscript.MS_ON
        self.layer.data= rasterpath


texts:
I maintain a Python package (Pyvisgraph) that builds a visibility graph given set of simple obstacle polygons. I am using it to calculate shortest paths between two points on a map. The obstacles are world shorelines obtained from here.

My problem is that the map is a "regular" projection (I'm a novice at GIS), so there is no visibility across the "seams" of the map. As an example, if I calculate the shortest path from San Fransisco to Shanghai, the path would not cross the Pacific, but go the long way through the Indian Ocean and Atlantic. There is no visibility from the US west coast to East China, because that is where the map boundary is.

I read up on gnomonic projection, but seems that would only work for one hemisphere at a time.

I am thinking that with any projection you will still have this "seam" problem. Is there any technique or projection I can use that would solve my problem?


texts: Update value in one table (feature class) by coping it from another table. The match between table is done on values from another field

texts:
I have two tables (actually both are feature classes in filegeodatabase).
They pretty much identical except. one is in lat/lon. second one is projected.
both have field line_length and objectid.
Second one has values in field line_length.
how do I copy these values (based on objectid field match) from second table into the fist one.

The snippet of my code:

Part first. Create dictionary from second dataset. Copy values from second dataset into dictionary. value that I will be using for matching is goint to be key, the value I need to copy to my first dataset are going to be values. This part works ok.

my_dict = dict()
fc = r'H:\NetworkAnalyst\TopologyVer4\firetrials.gdb\ft_mb_afterIteration_lambert'
fields = ['OBJECTID', 'line_length']
with arcpy.da.SearchCursor(fc, fields) as cursor:
    for row in cursor:
        a_id = row[0]
        b_len = row[1]
        my_dict[a_id] = b_len

print "dictionary done"

Second part. Update field "line_length" in the first dataset with values from second dataset, which are now in dictionary. Match records by OBJECTID field.

rows = arcpy.UpdateCursor(r'H:\NetworkAnalyst\TopologyVer4\firetrials.gdb\ft_mb_afterIteration_lambert')
for key, value in my_dict:
    match_value = key
    update_value = value
    for row in rows:
        print match_value, update_value
        if row.OBJECTID == match_value:
        rows.setValue("line_lenght", update_value)

I am getting an error for line
for key, value in my_dict:

Error is " TypeError: 'int' object is not iterable "

Any ideas?


texts: Why do XLS files cause Error: ".. is not valid or recognized data source"

texts:
In QGIS, I am trying to add some non-spatial data (in .xls file) that I want to use to join to some spatial data.

Adding the file, either from 'add vector' or dragging file onto project I get the above error.

Running QGIS v.1.7.4 on Windows.

Reading the documentation, this is supposed to work as of v1.7.

Wondering what I am doing wrong.

Also, tried to do a simliar thing with an MSSQL to a non-spatial sql table with the same error.


texts: For my undergraduate thesis I must do a spatial analysis over the OECD countries, but I cannot find a map for them, due to them being so spread apart. Funny thing, even my spatial analysis teacher couldn't do it, so I am wondering if there is any way you know about doing it?


texts: Problem with if-statement and whereClause format

texts:
I have a program which should count the number of a feature class or the length of a feature class depending on the geometry of the input feature class. This part is working.
There are two parameter-fields. One dropdown-field (optional) where you can choose a fieldname of the feature class and the other is a parameter field (optional) where you can choose a value depending on the fieldname which was chosen.
If I put in values in both optional parameter fields, the number or the length of the feature class should be count by using the inputvalues in a SQL-expression (whereClause).
If I am trying to do that, I get following error:
IndexError: tuple index out of range
in line:

whereClause = whereClause.format(arcpy.AddFieldDelimiters(fc_clipped, fieldname))

Second problem:
Why are both fields optional? Because if I don't put in any value to both optional parameters, the number of the feature class should be count (without a whereClause) with GetCountManagement.
But the code is still running through the else-statement and I get an error because of leaving the fieldname-parameter empty even though the parameter is optional.

So what I need is that the SQL-Expressions works with the input-parameters and
if the input-parameters are empty it should just use the GetCountManagement.

That is the full code:

    fc_clipped = os.path.join(out_gebietsscharfe_Objekte, objektname)+"_"+gebietsname #input feature
    fn = parameters[2].valueAsText #fieldname

    if arcpy.Describe(fc_clipped).shapeType == "Point" or "Polygon":
        if parameters[2].valueAsText == "#":
            COUNTER = arcpy.GetCount_management(fc_clipped)
        else:
            COUNTER = 0
            whereClause = """{} = '{}'"""
            whereClause = whereClause.format(arcpy.AddFieldDelimiters(fc_clipped, fn))
            for row in arcpy.da.SearchCursor(fc_clipped, fn, whereClause):
                COUNTER+=1
    else:
        fn = "Shape_Length"
        COUNTER = 0
        whereClause = """{} = '{}'"""
        whereClause = whereClause.format(arcpy.AddFieldDelimiters(fc_clipped, fn))
        for row in arcpy.da.SearchCursor(fc_clipped, fn, whereClause):
            COUNTER+=row[0]

    arcpy.AddMessage(COUNTER)


texts:
I'm using Leaflet for a project where I track the path of an entity along the map in realtime. I made a buffer with the last 1000 points but it is visible only through a big zooming cause the points are very near. So I'm thinking to develop an algorithm for building a smart path instead of keeping the last 1000 points but at the moment I'd like to know how many point can a polyline have in a normal PC. For example if I use 10K points, will the path continue to be efficient and fast? And with 100K? Is Leaflet optimized for polyline with many points?

A good solution could be to show only few points and load the path dinamycally when zooming (or maybe is it already done by Leaflet?).

I'm using Leaflet inside ReactJS.


texts: Error finishing importing Nominatim: idx_word_word_id?

texts:
After three weeks of hard work, I was one step from the end of importing planet.osm, but it crashed!

This is the first log error:

user@maps:~/Nominatim$ tail setup.log
  Done 37133684 in 148141 @ 250.664459 per second - Rank 30 ETA (seconds): -2.309861
  Done 37133685 in 148142 @ 250.662766 per second - ETA (seconds): -2.313866
  Done 37134058 in 151624 @ 244.908829 per second - ETA (seconds): -3.891244
  Done 37134058 in 151624 @ 244.908829 per second - FINISHED

Search indices
CREATE INDEX
connection to server was lost
ERROR: pgsql returned with error code (2)
pgsql returned with error code (2)

Then I restarted it with the comand written in the wiki, and then I got this:
    user@maps:~/Nominatim$ tail setup1.log
    Starting rank 29
      Done 0 in 0 @ 0.000000 per second - FINISHED

Starting rank 30
  Done 0 in 0 @ 0.000000 per second - FINISHED

Search indices
ERROR:  relation "idx_word_word_id" already exists
ERROR: pgsql returned with error code (3)
pgsql returned with error code (3)

and I obtain it everytime I use that command to resume the indexing. Don't tell me I do have to re-do everything from the beginning!

I use a Ubuntu 14.04.1 with the last PostGIS, PostgreeSQL, etc...

Can you help me?


texts:
I'm trying to use MODIS to download some spatial data, but I'm not able to do it. When I try to download and process MODIS data with MODIS R package, I get this strange strsplit error.
 b1 <- runGdal(product = "M.D13A2", begin = Sys.Date() - 30, tileH = 18:19, tileV = 4)

Returned:
########################
outProj          =  asIn
pixelSize        =  asIn
resamplingType   =  near
Output directory =  C:/*****************************/MOD13A2.061_20210525100758  (no 'job' name specified, generated (date/time based))
########################
Local structure is up-to-date. Using offline information!
Hm, I have to search for the file. Next time provide the full path and I'll be very fast!
Error in strsplit(x, ":") : non-character argument

I don't know what is happening here.


texts: Rename field in SHP file (DBF) with ArcCatalog

texts:
I want to rename some fileds on a SHP file, (the fields on the DBF file), using ArcCatalog (ArcGIS Desktop v10.2.2).

From ArcCatalog, right click on a SHP file -> Properties... -> Fields.
Then the program let you rename the field and change the data type... apparently.

When I hit Enter or click Apply, an error appears: Failed to alter the name of the field.

The file is not in use, nor read only, and the name of the field is under the 10 character limit.

Ironically, I can do it with gvSIG!


texts: Clipping raster inside a polygon (making hole)

texts:
Clipping a raster using a polygon with GDAL is easy - that is if you want to discard the data outside the polygon and keep what's inside it.
(Using gdalwarp -cutline)

I would like to discard what's inside the polygon and keep the outside (e.g. make a hole in my raster).

Does anybody know how to do this?


texts: GeoServer v2.9.1 re-sampling errors with GetCoverage request

texts:
App

I'm using Openlayers3 and GeoServer 2.9.1 to create an application allowing users to draw an aoi then download it as raw data subsets (geotiff) of a Canada wide data set (also geotiff). I'm using a WCS 1.0.0 GetCoverage request where bbox and Grid Dimensions are dynamically entered.

Issue

When validating the subsets, a simple subtraction of the downloaded data from the original dataset reveals pixel value changes (see image). How can I retrieve raw data subsets without re-sampling errors? It is critical that the pixel values be identical to the original.

Making a GetCoverage request from GeoServer's Demo Request Builder using the same bbox coords as my script gives me the same errors.

Pixel value difference from rasterCalc

Thoughts?


texts: Is there a tool to cut a polygon with the outline of a raster layer?

texts:
I would like my polygon layer to be cut into two where the end of the raster layer meets the polygon, is it possible? I've tried to polygonize my raster layer, but that doesn't seem to work. Any recommendations?
I'm working in QGIS 2.18.28.


texts: Mapserver 6.2 concatenate attributes in FILTER

texts:
Back in version 5.0 I had:

FILTER ("<tcode:expression/>" eq '[STATE][COUNTY]')
final version = FILTER ("11000" eq '[STATE][COUNTY]')

so I would filter out the features with state = 11 and county = 000

Now, in 6.2, it does not work anymore. Does anyone know how to make it work with concatenated attributes?

I know that it would work like this:

FILTER ("11" eq '[STATE]' and "000" eq '[COUNTY]')

but it would really help me to find a way for the concatenation to work like:

FILTER ("11000" eq ('[STATE]'.'[COUNTY]')) -no error, but not display anything

Been browsing the documentation for a while but I just can't make this work.


texts: Proper font size for each paper size for reader to be able to read names?

texts:
I am working on a project that I need to figure it out how or what font size for each paper size. There are 10 of them and two of them are in the same paper size. There is data frame I have set aside to display the crop name as a Index and need to display their names to idenfity them.

60 x 34
56 x 52
43 x 50
62 x 66
56 x 74

I am trying to put down all the crop names  for each polygon, for example here, .

I understand that it will be very diffcult to fill all the text names for each polygon for each paper size.  This one is the paper size 56 x 74.

I am curious if the reader would be able to read each font size at that kind of paper size? Any suggestions for me to tweak or fix them?


texts: Specifying different speeds for different vehicles when using ArcGIS REST API Vehicle Routing Problem service?

texts:
I am using ArcGIS REST API Vehicle Routing Problem service. I have several different types of vehicles, such as cars and trucks. There is an optional parameter called "impedance" to specify the travel times for vehicles. But it is very general, and it cannot specify for each vehicle.

What should I do if I would like different vehicles with different speed factors?

For example, cars with 1.5 while trucks with 1.


texts: Selecting closest points from a given point on another table

texts:
I have two tables, one with "houses", and the other with "points". I know how to get the distance individually, like this:

gisdb=# SELECT round(ST_Distance_Sphere(
(SELECT coordinates FROM houses WHERE id = 1),
(SELECT coordinates FROM points WHERE id = 2)))
as Distance
;
 distance
----------
      174
(1 row)

The thing is, if on the table "points" I have about 50 points, how can I make it return the distance of the closest one?

I've seen the function ST_ClosestPoint but I don't know how to apply it to my query.


texts: Lossless import of full OpenStreetMap planet dump into PostGIS

texts:
I would like to import a full OpenStreetMap planet dump as PBF into a PostGIS tables.

With tools like osm2pgsql or imposm3, configuration files are neccessary. Unfortunately, it is not possible to taking into account all possible attributes.

Is there a tool or way?


texts: Finding highest value within a shapefile or feature class using ArcGIS Desktop?

texts:
I want to find the highest value within a specific area of a shapefile feature class but am unsure about how to do this.

From the attatched image, I want to find the pixel within the outlined region with the highest value for fishing intensity (where the darker the rectangle, the higher the fishing intensity).

The background is a series of averaged point pixels rather than a raster layer


texts: Images in pop-ups on webmaps ArcGIS online

texts:
I am relatively new to webmaps - I am trying to add an image to pop-ups on web maps in ArcGIS and I am having problems getting the image to appear - all I can see is a little image icon with a tear through it as if it is having issues loading.

I have configured the pop-up media as they have listed in the instructions (http://doc.arcgis.com/en/arcgis-online/create-maps/configure-pop-ups.htm#ESRI_SECTION1_B16A9E11A48D48CBAFA9DB964FCEB805) and have replicated the images to a smaller size (thumbnails) around 200x150 pixels so they shouldn't be too big to load. I have also tried loading the URLs from Google drive.

The links actually work if you click on it and request to open the image in a new window - but i want the image to appear in the pop-up as well and it just wont appear - does anyone have any advice or solutions for this?

I read on a troubleshooting link that it may be because the url is a https rather than a http, would this be it?


texts: After clipping raster in ArcMap, file appears black in any other applications?

texts:
I have a tiff raster dataset, using arcgis I clip the data (to a feature class perimeter), the clipped file displays fine within Arcgis. However, if I try to view the clipped image in any external application the whole file appears completely black.

I have opened the clipped image in photoshop, gimp, windows viewer and paint and it always appears completely black.

I managed to get the image to appear by declaring the value of 255 (white) as the transparent colour, but now all the white detail has become transparent.

I am using ArcGIS Desktop 10.2.


texts: Is it possible to update column in QGIS using SQL statement?

texts:
In the current LTR of QGIS - 3.14 is it possible to directly update a column using SQL?
I understand in previous questions the advice was to use OGRINFO - Adding a column and update values using SQL in QGIS
Is there a way to run something like
Column_to_update =
{Select max(End_Chainage)
from DSC_Roads
group by RoadMntnc}

Currently the only way I can think of doing this is to run it (as below) in DB Manager, join by a key field and then update the column. Is there a single step to do the same?


texts:
I have a DXF file of a building from AutoCAD. I saved this DXF file as a SHP file in QGIS. But this SHP file is converted from a plain DXF file and is not referenced/projected, it is drawn from coordinates 0,0 in AutoCAD in meters.
And the drawing has meters, WGS84 has decimal degrees.

When I set Layer CRS as WGS84 in QGIS, the coordinates don't change.

How can I project this SHP file as WGS84 in QGIS?


texts: Get a Confusion matrix in QGIS from two vectors

texts:
How and with what QGIS tools can I get a confusion matrix from a vector polygon layer as a verifiable classification and a vector point sample layer as a reference classification?

That is, if the point with the class of "Rubber" hits the polygon c is a class "Rubber", the classification for this point is correct, if the point does not fall in any polygon, that point is not needed, if the classification of points and the polygon in which it are not the same, the classification for this point is wrong.

These are the vectors.


texts: Rules based rendering with spatial query

texts:
I would like to style polyline segments differently whether they are inside, on the boundary or outside a polygon. So based on the image above the red and green line segments inside the boundary would be styled differently

Some preliminary searches have shown that rules based rendering might be suitable but i'm not sure if I can use rules based on spatial queries like select by location.

My long way round would be to copy the layer, select by location, apply some filters and then style selected segments how I want them.

What I'm looking for is a way to automatically style those segments based on location so any line that goes inside that boundary get's styled in a standardized way.
 Currently using 2.14.2 Essen


texts: What is the standard method of finding satellite positions in between the 15 minute intervals in an SP3 file?

texts:
SP3 files provide satellite positions for every 15 minutes. How would you find the satellite position at a time between two logged positions, i.e. 7 minutes after the last logged position?

I can think of several ways to do this:
1. Just use the position at the closest interval.
2. Linear interpolation between the intervals on either side.
3. passing a polynomial through several consecutive points and using the polynomial to estimate any satellite position within the time frame.

I don't think the overall effect on the accuracy of the estimates would vary much between all three of these (but I don't know this as a fact). Is this the case?
Is there a "standard" way for dealing with it?


texts: How can we get a specific attribute of random position?

texts:
I am working on S57 chart. I want to apply *star shortest path algorithm on the S57 chart. Each node (position) should be defined whether is available for navigate or not. So at least each node should be on the sea.(whether on *SEAARE layer or not).

*S57 have different layer such as SEAARE which define sea water. And each layer  has OBJL attribute (SEAARE OBJL =119). This is specific integer number to define the layer.

  Pink area:Island - Blue area: sea. I want to detect by using API any random position on chart is navigatable(OBJL=119,blue) or obstacle(LNDARE OBJL=71,land)

I think that I can realize position on the sea using OBJL attribute but I don't know how I get these information from s57 for any position.

In short, I use QGIS and Python language for programming. How can I detect any attribute (like OBJL) for any random position (lat/long) on chart by using pyqgis or GDAL/OGR API?

If you have any idea of other file format, I can convert s57 to shapefile or PostgreSQL. Also I am quite new to GIS.


texts: Efficiently converting/storing WGS84 Lat/Long as Integers

texts:
I'm a GIS newbie and I'm trying to figure out how I might store WGS84 coordinates as an integer efficiently. I could simply multiply it (10e7 for example) but I feel like this must be a problem that's been solved before. It seems, though, that there isn't any standard that uses integers.
I'm using datasets that require integers for co-ordinates so need to create my own integer co-ordinate system, and I'm struggling to "stand on the shoulders of giants" with this particular problem.
Is there some standard to do this?
It looks like google just uses 10e7 for their tracking information for example.
Is that the most efficient way to do this?


texts: Preventing shapefile lines becoming rasterized when exported as PDF by ArcGIS for Desktop?

texts:
In ArcGIS for Desktop some Line feature shapefile becomes raster when exported as pdf.

In exported PDF the line features layer is not present as it becomes raster.

The line features are mixed with satellite image and not exported as different layer in PDF.


texts: How to find the size of the ArcGIS map when the map takes up less space than the containing div (zoom=0)

texts: I've been digging through the ArcGIS map trying to figure out how to find the size of the map. The map object has width / height but those are the containing div and not the width / height of the current map. Is there a way to find / calculate the width / height of the map?


texts: Applying Stretch to Multiple Rasters in ArcGIS Desktop?

texts: In ArcGIS Desktop 10.0, how can I apply a stretch, which I would normally do by selecting 2 standard deviations in the "Symbology" tab of Layer Properties, to 200 rasters simultaneously?


texts: Adding a new Shapefile to GeoServer through REST API

texts:
I'm trying to interact with GeoServer's REST API using PHP's cURL and I'm having an issue when uploading a new Shapefile...

This is my uploadShapefile method:

public function uploadShapefile($zip, $ws, $ds) {

    if (!$ws || !$ds) {
        return false;
    }

    $fd = fopen($zip, 'rb');
    if (!$fd) {
        return false;
    }

    $headers = ['Content-Type: application/zip'];

    $settings = [
        CURLOPT_HTTPHEADER     => $headers,
        CURLOPT_USERPWD        => $this->passwordString,
        CURLOPT_PUT            => true,
        CURLOPT_UPLOAD         => true,
        CURLOPT_REFERER        => true,
        CURLOPT_BINARYTRANSFER => true,
        CURLOPT_INFILE         => $fd,
        CURLOPT_INFILESIZE     => filesize($zip),
    ];

    $endpoint = "rest/workspaces/$ws/datastores/$ds/file.shp";
    $response = $this->call($endpoint, $settings);

    fclose($fd);
    return $response;

}

My call method just uses cURL and returns the HTTP Response code

private function call($endpoint, $settings) {

    // $this->service is 'http://localhost:8080/geoserver/'
    $ch = curl_init($this->service . $endpoint);

    foreach ($settings as $setting => $value) {
        curl_setopt($ch, $setting, $value);
    }

    // Execute
    $buffer = curl_exec($ch);

    // Result
    $info = curl_getinfo($ch);
    return $info['http_code'];
}

When I call the uploadShapefile method it creates the datastorage in the specified workspace, BUT its Shapefile Location under Connection Parameters is empty and it also doesn't create the new layer for the feature type.

I've checked my $fd and it can read the zip file. Am I missing something?

Edit:

If I use CLI curl the behaviour is the same (datastorage created but no layer) and I get:

*   Trying ::1...
* Connected to localhost (::1) port 8080 (#0)
* Server auth using Basic with user 'admin'
> PUT /geoserver/rest/workspaces/myworkspace/datastores/streets/file.shp HTTP/1.1
> Host: localhost:8080
> Authorization: Basic YWRtaW46Z2Vvc2VydmVy
> User-Agent: curl/7.49.1
> Accept: */*
> Content-type: application/zip
> Content-Length: 215432
> Expect: 100-continue
>
< HTTP/1.1 100 Continue
* We are completely uploaded and fine
< HTTP/1.1 202 Accepted
< Date: Mon, 14 Nov 2016 17:26:40 GMT
< Date: Mon, 14 Nov 2016 17:26:41 GMT
< Server: Noelios-Restlet-Engine/1.0..8
< Transfer-Encoding: chunked
<
* Connection #0 to host localhost left intact


texts: Trying to import Shapefile to SQL through Command

texts:
I am a newbie and I'm trying to import a shapefile into SQL Express 2014 via the OGR2OGR tool and I am facing errors.

I installed SQL Express 2014 following this article.

http://kb.printaudit.com/display/PA6/SQL+Server+2014+Express+Installation

Afterward, I followed the below-mentioned article for Shapefile Import.

https://alastaira.wordpress.com/ogr2ogr-patterns-for-sql-server/

Then I installed QGIS version 2.10. This is the command I entered:

ogr2ogr -f MSSQLSpatial "MSSQL:server=JDOE\SQLEXPRESS;database=MySpatial;trusted_connection=yes" C:\Users\jdoe\Desktop\ShapeFile\MyShapeFile.shp -nln "Spatial" -progress

This is the error I am getting:

ERROR 1: INSERT command for new feature failed. [Microsoft][ODBC SQL Server Driver][SQL Server]Arithmetic overflow error converting numeric to data type numeric.
ERROR 1: Unable to write feature 0 from layer MyShapeFile.
ERROR 1: Terminating translation prematurely after failed translation of layer Ontario_FED_2015 (use -skipfailures to skip errors)

Please suggest what am I doing wrong here.


texts:
I am trying to address some latency issues accessing US based servers from Asia

We are using GeoWebCache in front of two GeoServer instances, based in the USA.

One thought that occurred to me was to deploy additional GeoWebCache Servers, but instead of having them hit a GeoServer instance, locally or remotely, have them front another GeoWebCache instance (Cascaded) located in the USA.

The thought behind this is that the tiles will already be cached on the "Master" GeoWebCache which simply serves the requested tiles to the GeoWebCache instance requesting it in Asia. This should eliminate the need to have additional Geoserver Instances to service Tile requests

Has anyone done this before or have any thoughts as to whether it will work as I expect?


texts: How to make dynamic labeling in a project for new features on QGIS?

texts:
So let's say I have an integer id column that users regularly add features to. mytable1 has 500 features and growing.

When doing a categorized view on these, I'll get 500 different colors for all these features. Now problem is that when a new feature gets created, it does not get styled in any way. I can introduce a new style with an else but then every new feature will be the same color

Now is there a way that every feature with an id ending in 0 to be red (essentially, id 0, id 10, id 110,...) and every feature with id ending in 1 to be blue (id 1, 11, 31, 61,...).

Bonus if I was allowed to actually have 30 colors instead of just 10.

So ids 1, 31, 91, 301 will have same color

And ids of 29, 59 , 89 will have same color

Is there a way of doing this in QGIS?


texts: How to normalize Landcover criterian Map to 0-255 scale in Weighted Linear Combination?

texts:
I am interested to normalize Landcover raster image having 6 classess using fuzzy membership function (sigmoidal monotonically increasing) in Weighed Linear Combination (WLC) in the backdrop of MCE approach.

The aim is to produce Landslilde susceptibility map. The weights will be derived using pairwise comparison of the attributes to get their relative importance over one another.

Normalization/standardization of input criterian maps is a prerequisite in WLC. I have already ranked my classes from 1-6. Let me know how could I assign inflection points to my landcover classes on 0-255 normalization scale to produce my fuzzy-landcover map?


texts: Seeking Detection Tutorial for Google Earth Engine

texts:
Does anyone know a good tutorial to learn how to detect objects using Google Earth Engine?
My goal is to find a smart way to detect objects (these objects are plantations in Ivory Coast) using Landsat data.
Here is one example of the plantation I am trying to detect : A Plantain Plantation near Anyama 5.464013217554387, -3.8705557721123736


texts: UID for set flag and barrier commands in arcobjects .net

texts:
I have found below UID to add a flag and barrier.

Flag and Barrier Tools : {C86085FD-D8D9-11D2-9F4F-00C04F6BDF0E}

But nothing is executed when using this UID. I would need to know the UIDs for junction flag, junction barrier, edge flag and edge barrier to invoke them from my application developing using arcobjects .net.


texts: Get ESRI JSON response locally from Arc document?

texts:
I'm trying to use the ArcGIS Flex SDK for a mobile app that will sometimes be used in environments with no connectivity. I've figured out that layers can be brought in locally by saving the JSON response given by the server and simply using the local file as the layer URL.

Is there a toolbox, plugin, or .py script that will take existing features or shapefiles and convert them to the ESRI JSON format that the server uses? I've tried esri2open, but that only does GeoJSON and a generic JSON object, and I'd prefer to have something that automatically gives the exact same JSON structure as the REST services would when loading an online map. Ideally, it would be something that could run from the Arc desktop without the additional steps of publishing a service and downloading the response. Does something like this exist? I've noticed that 10.2 may have that capability, but for now I'm stuck with 10.1.


texts:
I have a CSV file with mainly UK data - postcodes that I have then referenced to there Latitude and Longitude. I have 5 rows that are outside the United Kingdom, so I have looked up the Latitude and longitude of these points and also added these to the dataset.

I am trying to plot these points on a world map - the aim is to eventually create a webmap using qgis2web. I have got a shapefile I have downloaded from Eurostat:
1:10 Million    World   Point / Line / Polygon  Shapefile   2014    ETRS89
When I add my delimited text layer in QGIS using the latitude and longitude values the points are way off where they should be even when choosing ETRS89 as the coordinate system.

The majority of the points should be sitting above London in the UK, but are actually off the coast of Africa!

I have tried using different coordinate systems and basemaps, but I cannot get any of them to line up. I know it is probably something really simple, but it is baffling me currently.


texts: How to use GROUP BY with arcpy.da.SearchCursor in file geodatabase?

texts:
fc = "C:\Users\AuxS3-4\Documents\ArcGIS/Default.gdb\Malla"
    fields = ('OBJECTID','Sentido','Carriles','idMalla')
    sqlclause=(None, "GROUP BY Sentido")
cursor = arcpy.da.SearchCursor(fc, fields,sql_clause=sqlclause)

not working

  RuntimeError: An invalid SQL statement was used. [SELECT
  OBJECTID, Sentido, idMalla FROM Malla GROUP BY Sentido

arcgis 10.1!!!


texts: Is it possible (and how) to mask NoData values using Scikit-learn PCA?

texts:
I have to carry out a Principal Component Analysis (sklearn.decomposition.PCA) on a dataset composed of several images. Some of these images have very large gaps due to the removal of clouds.
Is there the possibility to apply a mask to all the NoData values for each layer so as the algorithm does not take them into consideration, and consider only the valid values of other underlying/overlying layers (which do not have NoData in that position)? So to have final PCs without gaps.
I'm not sure if my request is clear, I try to give an example:
I have a dataset consisting of 4 layers. In the pixel in position (x = 27; y = 32) the different 4 layers have:
Layer-1 = value;
Layer-2 = value;
Layer-3 = NoData;
Layer-4 = value.
I need that the alghoritm takes into consideration the Layers-1,2,4 and discard the Layer-3 to build the PCs values of only that pixel.