texts: Use PYQGIS to view MS SQL Server temporary table in QGIS

texts:
I'm working on a QGIS plugin to query a non spatial database held in MS SQL Server and to display the results of the query as a non spatial table in QGIS. I've tried the following:

create global temporary table (I've simplified the select query in this example):

sql=("SET NOCOUNT ON;"
        " IF OBJECT_ID('tempdb..##t1') IS NOT NULL DROP TABLE ##t1;"
        " USE MSSQLDatabase;"
        " Select * into [dbo].##t1 from MyTable;")
 conn = pyodbc.connect(con_str, autocommit=True)
 c = conn.cursor()
 c.execute(sql)

display the table in QGIS:

uri = QgsDataSourceURI()
uri.setConnection(".\SQLEXPRESS", "", "tempdb", "", "")
uri.setDataSource("dbo", "##t1","")
vlayer = QgsVectorLayer(uri.uri(), "t1", "mssql")

The code to create the global temporary table works - I can query the temporary table in another app, but my code to display the table in QGIS doesn't produce a valid layer. The QgsDataSourceURI code works if I point it to a non temporary table in the database, but I can't get it to work with a temporary table. Hope someone can spot what I'm doing wrong!


texts: Getting full path for feature classes in dataset using ArcPy?

texts:
I have a python script that I am writing to create a copy of a database and then clean it up (basically remove any empty items). Everything is working with the exception of feature classes that are embedded in a feature dataset. I have the portion of code that iterates through these items and is the source of my issue below:

arcpy.env.workspace = inputGDB
dataset_list= arcpy.ListDatasets()

for dataset in dataset_list:
    featureclass_list=arcpy.ListFeatureClasses("","",dataset)
    for feature in featureclass_list:
        inpath=arcpy.Describe(feature).catalogPath

When I do this, the "catalogPath" returned is not the full path. It only contains the path up to the gdb and then appends the feature class name to it (e.g. it returns D:\Folders\Database.gdb\feature_name not D:\Folders\Database.gdb\Dataset\feature_name). The dataset which this feature was originally in is not returned with the above code. Now if I do:

item=arcpy.Describe(r"D:\Folders\Database.gdb\Dataset\feature_name")
print item.catalogPath

I get the correct path for the data in question. So the question I have is how do I get the full data path for feature classes embedded in a dataset with some option similar the the first snippet of code?


texts: Save button and TOC layers missing from Portal Web Map (ArcGIS Enterprise)

texts:
I have a few web maps in portal that correctly show data without any errors but there are no layers shown in the TOC (Contents) and the Save button is missing.

This sometimes happens when I update a service, more specifically, delete a service and publish a service of the same name (using Pro "Map Image Service" (cached) since these cannot be simply overwritten as vector services can)

If I remove the layer for that particular service from the web map prior to deleting it, then re-creating it, then add to the map, everything works fine. If I don't do that the TOC (Contents) no longer shows layers and the save button is gone. If I search for layers in this (now broken) web map, it still works and if I add a new item, the TOC layers mysteriously appear! However, since the Save button is gone, I can't repair and update my web map by simply saving it.

When I use the the ArcGIS Online Assistant to repair any broken links it doesn't really do anything since the links are valid and are actually identical to what they were before. Furthermore, the new "updated" service correctly appears in the web map, just the TOC has no layers and the Save button is missing.

I should add that I am the admin and have all admin privileges. I looked online for some answers but only found one post on GeoNet but none of the suggestions there offer any solutions or workaround.

How Can I bring back the Save button?

Example of a "broken" web map (TOC layers missing and no save button):

Example of a correct web map (TOC layers and button shown):


texts: QGIS v3 ImportError: No module named builtins

texts:
I am trying to automate some tasks in QGIS, and I was initially going to use the model builder but I thought I would try some Python instead, using PyQGIS.

I am using this as a quick start guide as I am relevately new to Python.
https://www.qgistutorials.com/en/docs/running_qgis_jobs.html

Anyway I got the "ImportError: No module named builtins" when trying to run a batch file.
Having looked on this forum I seen others have the problem (which was a relief) and I see it is a problem with the new version of QGIS, and the solution seemed to be running a batch file in QGIS 3.0/bin/python-qgis.bat which sets up the paths correctly, however I have run this then tried re-running my batch file, and the error still remains.

Not sure if I am just doing something wrong, I am quite new to coding so sorry if it is something obvious.


texts: How can I get SOESupport to handle Polymorphism?

texts:
I've got some business objects below that I can serialize and deserialize successfully with Json.NET.  I don't see a simple way to do this with Esri's SOESupport however.  I'm happy with Json.NET, but I'd like to leverage Esri's IGeometry conversion routines.  Is there some way to get SOESupport to work for this?

Update:

Since it doesn't look like Esri can do this, I'm trying to use Json.NET. See my related issue on stackoverflow:

How can I deserialize with TypeNameHandling.Objects in Json.NET Silverlight?

public class Animal
{
    public string Name { get; set; }
}
public class Dog: Animal
{
    public bool LikesBones { get; set; }
}
public class Cat: Animal
{
    public bool LikesMice { get; set; }
}
public class Zoo
{
    public int ID { get; set; }
    private List<Animal> m_Animals = new List<Animal>();
    public List<Animal> Animals { get { return m_Animals; } set { m_Animals = value; } }
}

Is there some way to get Esri's JsonObject to support this sort of test?

public static void Test1()
{
    // use Json.NET
    Zoo z1 = new Zoo() { ID = 1 };
    z1.Animals.Add(new Dog() { Name = "Fido", LikesBones = true });
    z1.Animals.Add(new Cat() { Name = "Felix", LikesMice = false });
    string s1 = Util.Serialize(z1);

    Debug.Print(s1);
    var z2 = Util.Deserialize<Zoo>(s1);
    string s2 = Util.Serialize(z2);
    Debug.Print(s2);
    if (s1 != s2)
        Debug.Print("failed");
    else
        Debug.Print("ok");
}

I can apply settings in Json.NET that force the typename to be serialized, allowing polymorphism for the list of Animals.  How can I implement these same functions with SOESupport?

public static string Serialize(object o)
{
    JsonSerializerSettings settings = new JsonSerializerSettings();
    settings.TypeNameHandling = TypeNameHandling.Objects;
    return JsonConvert.SerializeObject(o, Formatting.Indented, settings);
}
public static T Deserialize<T>(string s)
{
    JsonSerializerSettings settings = new JsonSerializerSettings();
    settings.TypeNameHandling = TypeNameHandling.Objects;
    return JsonConvert.DeserializeObject<T>(s,settings);
}


texts: How to make Map Book of pages with same area shown but different layers turned on/off?

texts:
In ESRI a Map Book is best suited for a series of features that each get their own page. It allows you to pick a feature and generate a Map Book Page for each feature including any additional information you wish to add to each page via dynamic text, charts, etc.

Is there an alternative for say a series of pages showing the EXACT same area but with different layers turned on and off? We have a Capital Improvements document that I have to create showing the City Boundaries on every page. One page has all police and fire stations listed, another shows proposed water systems, another shows expressways and highways and work being done near or on them... etc., ad nauseum.

Is there some way to set this up so that I DON'T have to created a separate MXD for each Map Page of this document?


texts: Teaching GIS, what to do and what not to do?

texts:
Because I will (in the near future) be a teacher of GIS. I need some advice form people who have previously done this. Any clues - what to do and what not to do ?
I think it will be very helpful for my preparation.


texts: How do I change the colors of my layer file?

texts:
I have a ILayer which I load from a .lyr file which I created from a shapefile.  How can I change the colors of the elements?  I want to make the roads red instead of green for example.

I need to do this programatically with ArcEngine 10 and C# VS2010.


texts: Is there a way to control which ArcMap AddIn is loaded first?

texts:
I have a shared DLL that use across multiple ArcMap AddIns.  It rarely changes, but recently I've had to modify it and I started getting MissingMethodExceptions whenever I made calls to classes in that shared DLL.

After some investigation, I realized that there are other ArcMap AddIns that carried an older version of the shared DLL that were being loaded first. And if I clear out all of the other AddIns and blow away the assembly cache, I run by AddIn with new Shared DLL with no issues.

Is there a specific order of AddIn loading from the assembly cache? In my case, it appears to be oldest first. Is there a way to alter the order of loading?

Related Links:
http://forums.arcgis.com/threads/29870-Add-ins-that-share-same-custom-library


texts: Fill/interpolate no data values on raster

texts:
I have a large stack of rasters for snow cover every day for the past fifteen years. Each pixel can have one of three values: snow, no snow, or No Data due to cloud cover.

I want to fill those no data values with data from the following days, i.e. If 12/1/2002 has a blank spot due to clouds, fill that spot with the data from the 12/2/2002 raster. I do not want to alter the other pixels at all- just the blank/No Data pixels.

I'm using ArcMap 10.1

Is this an approach which will fill the No Data values?


texts: Couldn't load QGIS utils; Python support disabled; No module named future

texts:
Unexpectedly Qgis won't open today. Error code is:

Couldn't load QGIS utils.
Python support will be disabled.

Traceback (most recent call last):
  File "", line 1, in
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 20, in
    from future import standard_library
ImportError: No module named future

Python version:
2.7.10 (default, Oct 14 2015, 16:09:02)
[GCC 5.2.1 20151010]

QGIS version: 2.15.0-Master 'Master', exported

Python path:
['/usr/share/qgis/python', u'/home/simon/.qgis2/python', u'/home/simon/.qgis2/python/plugins', '/usr/share/qgis/python/plugins', '/usr/lib/python2.7', '/usr/lib/python2.7/plat-x86_64-linux-gnu', '/usr/lib/python2.7/lib-tk', '/usr/lib/python2.7/lib-old', '/usr/lib/python2.7/lib-dynload', '/usr/local/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages/PILcompat', '/usr/lib/python2.7/dist-packages/gtk-2.0', '/usr/lib/python2.7/dist-packages/ubuntu-sso-client']

I've done some digging around similar errors to no avail. utils.py imports a module from __future__ successfully then standard_library from future unsuccessfully. This suggests that future/standard_library should work though. I've just reinstalled (xubuntu, debian nightly ppa) and still no joy. Qgis has been working fine until now, although I've not opened it for about a month or so. Searching for future in /usr/lib doesn't give any other exact matches: futures.py in /python3.4/asyncio, fix.future.py in python3.4/lib2to3/fixes, the same file in python2.7 folder...

2.7 isn't new, so 2.7's future has assumedly never had standard_library in it, right? What am I missing here? Could an independent package have decided I no longer need a more recent version of python, I okayed it for autoremove, and it turns out Qgis needs that version of python's future?

OS is xubuntu 15.10


texts: Measuring the distance between lines and points in QGIS

texts:
Using QGIS - I have a vector of fault lines.  I have point data in the form of a csv which I have also loaded into QGIS.  I need to know how far each point is from the nearest fault line then enter this distance in a new column in the csv file.
Any ideas?


texts: QGIS Polygonize (GDAL) problem - Null value in attribute

texts:
Is it possible to use Polygonize in QGIS (GDAL) so that a line shapefile is converted to a polygon shapefile with an attribute(s) in the line file carried into the polygon file.  I am trying to build an elevation polygon dataset  - I downloaded data from EU-DEM, clipped and merged a few files into a new raster.  Then created contours from the new raster.  Now I want to create polygons - however the elevation value is not retained in the output polygon shapefile.
I have said yes to the option Keep Table Structure of Line Layer.
THe struture of the polygons are all fine but the value in the elevation attribute is NULL.
One way around of this is to do it for each value of elevation 100, 200....up to the max.  Then add an attribute and set it to the correct value - however it would be quicker if Polygonize could do it.  Am I missing something here?  Is this a bug in Polygonize or was it just not designed to do this? Anyone got any ideas?  I am using QGIS 2.4.


texts: I am trying to convert a raster to polygon in ArcGIS 10.2.2 but when I run the tool, it says complete but no polygon file is created and the tool window pops up again. I have tried converting the raster to integer type, but still no luck! I am using a raster projected from Mollweide projection to GCS WGS 1984 and both the input raster and output location is in a geodatabase.


texts: Query all Object within n kilometer Radius with Hibernate Spatial?

texts:
I use hibernate spatial attach a geolocation to a car. My card domain class looks like this:

import com.vividsolutions.jts.geom.Point
class Card {
  String name
  Point location
}

My program is in Grails so the samples I gave are in Groovy. I found a similar post here which does not really answer the most important question on how to specify the radius correctly to set n kilimeter for the radius.

Here is how I compute a circle Geometry:

  private static Geometry createCircle(double x, double y, final double RADIUS) {
    GeometricShapeFactory shapeFactory = new GeometricShapeFactory();
    shapeFactory.setNumPoints(1000);
    shapeFactory.setCentre(new Coordinate(x, y))
    shapeFactory.setSize( (RADIUS * 2)/88.1)

    return shapeFactory.createCircle().getBoundary()
  }

The size of the circle is divided by 88.1 the is just a dirty fix to get an approximate dimension but it is still wrong.

My query is done like this:

double radius = 40
Geometry filter = createCircle(car.location.x, car.location.y, radius)
Session session = sessionFactory.currentSession
Query q = session.createQuery("select c from Car c where within(c.location, ?) = true")
q.setParameter(0, filter, GeometryUserType.TYPE)
q.list()

This works not very accurate. Some of the cars which should be outside circle are returned from this query.

Here is an example. My center is Hamburg and the radius is 40km. I made a google maps visualization.

Here is when I set radius = 40:

You can see that at the top left one car which has a location outside of the circle is still drawn. This should not be the case. It appears to me that the circle I draw with google maps is not equal to the circle Geometry I draw in code for my query.

Here is when I set radius = 30:

You see that the cars at the bottom right disappear which is correct, but the car at the top left still remains in the query.

When I draw the circle I created with createCircle I get the following (using getCoordinates() to get the coordinates of the circle):

How can I query all cars within a radius of 40km? How do I get the correct coordinate system?


texts: Iterative executing select layer by attribute

texts:
I have a Shapefile. There  is a field "YEAR" (range 2000 - 2015) in the attribute table.

I want to select layer base on different year and save as a new layer and don't know how to add a loop into following model.

I try to add a for loop.

The error occurs.


texts: QGIS, What is the reason for wrong scale and wrong areas?

texts: I am currently working with QGIS and loaded openstreetmap as a basemap (web-> OpenLayers plugin), it is in WGS84/ Pseudo Mercator. For some reason the scale is not correct. For example a distance in my project results about two times as long as the same distance for example on google maps. Does anyone know what the reason could be?


texts: I'm using mapbox to show gravestone locations for a cemetery. Is there a way to have the list on the left be split into 5 items, with subcategories in each? IE: Last names, "A-D", "E-I", etc. and then the actual persons name in the subcategory list?


texts:
I have seen Downloading Subarea vs Clipping Images to AOI using Planet API? and I still do not feel satisfied with the answer so raise this again.

So what is really the difference between them? @bosth in their answer said

  An important distinction is that the Clip API will be decomissioned in
  favour of other means of downloading parts of images, but the
  documentation is still there for existing users to reference.
  Currently, it is recommended that new users employ the subarea method.

I think the difference is Clip to AOI clip right at the server side, at the API endpoint https://api.planet.com/compute/ops/clips/v1/ and you not only download small size image but also pay for small size footprint/area why with the Subarea you pay for the whole scene footprint/area of 220km2 even if you clip only 1km2 out and download. In my case, when I activate the scene, get the download link, pass it to VSICURL and use gdal.warp() Planet actually calculate the area of whole scene and accumulate in my Used Area. The only thing good with Subarea is that you download the small size image which is faster, low bandwidth and reduce your local storage.

And if Clip API will be decommissioned (or has already!) then the only way to clip and download is Subarea, which seems not to save your Km2 quota.

Does anyone know the correct answer/explanation for this case?


texts: How to edit a model that is saved in project?

texts:
I have a model saved in QGIS Project which i can find under "Project models". But how can I open and edit it?

I have tried dragging and dropping it to the empty canvas, and using open Icon in toolbar (only .model3 files can be opened there), double clicking (runs the model), right clicking etc.


texts: Water Utility Network Tool - Editing the Configuration Files

texts:
Recently I've obtained the Water Utility Network Tool developed by ESRI. (Found here: http://solutions.arcgis.com/utilities/water/help/network-editing/ )

In working with it, I've attempted to visualize jump overs/unders when a water main crosses over another without intersecting. Which this tool does exactly as stated, except the actual visualization of the jump is tiny and when zoomed out beyond say 1:15, you cannot see it anymore, yet for purposes of my edits, I would like to see it on the map visually in a smaller scale setting. (For sizing representation, the map will show a multi-city block area per map, whereas the current jump radius is only visible at the sub block level).

I have found the configuration.xml files within the tool via 7zip archive. The first file Config.xml at the start appears to only handle how the button looks and where to find the image but not control the radius. I then found the second file under the 'install' folder called loaded.water.config which has this line of code within it

<!--The radius of the semi circle created by the Jumper tools-->
<add key="CreateJumps_Distance" value="0.5" />

This is the default setting. I've tested changing the value to 2.0 and even up to 5.0 and then saving the file and having the archive update.

Problem is, nothing changes in ArcMap Desktop then, it remains as if .5 is in effect that the radius didn't change. I attempted to re-add the add-in after the archive update but states that its an invalid digital signature now and won't install saying "Installation Failed. Access is Denied. The specified add-in may be in use." even though I've closed the 7zip and ArcMap.

Am I not updating the proper area to make the radius of the jump circle larger? Missing a step on updating this code to work?


texts: Creating 2 symbology in the same layer on Mapinfo

texts: New to MapInfo, I would like to create a High voltage layer with two different attributes (11 and 6.6) representing by two different colours on MapInfo. How can I do that?


texts: Convert points(.shp file) to raster (.tif) in QGIS

texts:
I have some points in .shp vector file and I need to convert them in .tif file with 1km x 1km resolution. The problem is - there are some points in one cell during conversion, and some values are missed. How can I convert these point without missing? Is there any function that will sum values in cell(in QGIS)?
Any suggestion are acceptable.


texts: GeoWebCache not deleting expired/discarded tiles from _gwc_in_progress_deletes_ folder

texts:
Is it OK to 'force' delete folders/files from GeoWebCache's

_gwc_in_progress_deletes_

folder?

I have a folder (zoom level) with several subfolders each having tens of thousands of tiles that have not been deleted for about a month now (& with no detectable changes/no errors in geoserver error log).

GeoServer Version 2.12.2

GeoWebCache version 1.12.2, build 1.12.x/8e746f758aabae3ab832c30757de6fafa4c39d0a


texts: How to determine maxResolution and numZoomLevels in dynamic environment?

texts:
I want to determine resolution and number of zoom level dynamically instead of inputting the values in openlayers client. When the user move the position on map by dragging or panning control, I can track the updated Lat-lon and get the min and max lat-lon. Now, I would like to calculate the resolution and maxzoomLevel for newly calculated latlon position when each time user changes the map position.

Thanks


texts: How can I use use GRASS that was installed from an image in Python3 on MacOS?

texts:
I am running MacOS 10.15 and want to run GRASS 7.9 on my machine from a Python program. So I downloaded GRASS as an image here.

Now I want to import the GRASS functionality into a Python program.

I have tried going into the resources of the app in my applications folder and add the Python GRASSs folder to my path. I have installed grass-session and tried to set the GRASSBIN variable to the application path. But so far I have had no success.

I have installed grass-session version 0.5 with pip, set the GRASSBIN environment variable to /Applications/GRASS-7.9.app/Contents/MacOS/Grass and ran import grass_session, but I am getting the following error:

RuntimeError: GRASS GIS start script: /Applications/GRASS-7.9.app/Contents/MacOS/Grass --config path, return as GISBASE a directory (tab 1 of window id 16042) that do not exist.

and GRASS opens itself!

Does someone know how to use GRASS in Python?


texts: Checking if raster extent overlaps AOI in ArcPy?

texts:
I have have a directory that contains many of rasters. I also have a polygon shapefile with a single record that defines my area of interest (AOI). I would like to copy all rasters that fall within (or partially within) said AOI to a new directory)

Currently, my code is like this

#set up paths and get list of rasters
arcpy.env.workspace = r'\path\to\raster\dir'
rasters = arcpy.ListRasters()
AOI = r'\path\to\aoi.shp'

#get geometry of single polygon in shapefile
with arcpy.da.SearchCursor(AOI,['SHAPE@']) as cursor:
    for row in cursor: AOI_geom = row[0]

for raster in rasters:
    raster_extent = arcpy.Describe(raster).extent
    if raster_extent.overlaps(AOI_geom):
        #copy raster to new dir....

However, the extent.overlap() always returns false, even when I can see visually that there is an overlap


texts: Converting single point layer to buffered raster using QGIS?

texts:
I'm trying to use the proximity tool to create a distance raster from a single vector point. I found here: Calculating distance to points in QGIS and [d-3] Create Proximity Raster. However, I couldn't go past the first step which is rasterizing the vector point layer.

I couldn't figure out how to rasterize the single point and extend the raster to a specified distance, like a buffer, and set all the other values to nan or zero. I've tried the different rasterizing tools available in QGIS but I couldn't find anything that works.


texts: Running a python script that will find fields that are not attributed

texts: I have about 30 geodatabases with about 15 feature classes in each database.  Is there a way to have a script check the attribute tables to see if they're not attributed?  If it could then select the features that don't have attributes and export them to a gdb or shp.  I want an easy way to check the tables for nulls that I could run nightly.  Would a searchcursor work for this?  Is this possible?


texts: Creating point layer from linear layer with only ends of lines using ArcGIS Desktop?

texts:
I have ArcGIS 10 with ArcView licence and I need to make point layer from a linear layer, with only start and end points for each line. In ArcInfo I could use Feature Vertices To Points with "both_ends" as option but I don`t have such licence.

I found "Polyline Vertices to Points" script for creation point layer from linear layer with only ArcView licence, but it returns all vertices not only ends.

How can I delete all vertices from lines but ends or maybe other variant how can I make point layer from line layer with only ends of lines?


texts: How to get a subset of Mathew Hansen dataset in projected co-ordinate system

texts:
As we know from the metadata of the Global forest datasets created by Mathew Hansen, et. al., that can be found here, this data is derived from the landsat data. So the spatial resolution of this data should be 30 by 30 meter.But the problem is they are not- or maybe I am on the wrong boat.

This data comes in WGS84 projection system but when I reproject subset of this data into PCS(UTM46), I do not get the 30 by 30 resolution rather ~25m.

More dangerously, when I just extract by mask this dataset to my region of interest it gets displaced as the below visual.The mask featureclass is of same projection system i.e. WGS84 just like the Hansen dataset.

Now can anyone help me with what the correct/scientific way of getting raster dataset for specific area(subset, UTM46) in projected co-ordinate system from this Hansen dataset?


texts: Downloading Sentinel-1 SLC in individual bursts

texts:
I'm about to embark on a DInSAR time series project, for which I will be using Sentinel-1 single look complex (SLC) data.
Sentinel-1 SLC products are only available in Interferometric Wide-swath (IW) format, which consists of three subs-swaths, further subdivided in bursts (see figure 1 on page 3 of this document). The data volume for the complete time series I want to analyze would be around 8TB when downloaded with complete SLC images. Since I don't have the data storage capacity for these volumes and given that my research area is rather small, I would prefer to download individual bursts per SLC image.

Would it be possible to download individual bursts of an SLC image as opposed to the whole wideswath? If so, how can it be done with the Copernicus Scihub API?


texts: Display map tip in QGIS does not show latest image while hovering

texts:
In QGIS 3.14 project I have a polygon layer with fixed names (zones) and I generate png figures (with pyqgis script) with the name of the zone to show them while hovering with map tip (Display).
However, the popup figures are refreshed only in some polygons, when I replace these png plots and in some polygons, old plots are shown which were already replaced. Is there a way to force QGIS to show latest files?
Here's HTML I use to show the figures:

I tried to trigger repaint or refresh all layers with the following code with no result.
layer.triggerRepaint()
canvas.refreshAllLayers()

adding '_plot.png'%]?cache=none"/> also didn't help
Note: I currently found a workaround by adding new field where I add a random number and also add this number to figure name and use it to generate path, but it doesn't look like the best solution...


texts:
Im trying to perform a loop on the Select_Analysis tool in Arcgis. Code looks like this:

inputfile = open("locality.txt", 'r') # locality.txt is a txt file of locality names
loc_list = []
for line in input_file:
    loc_list.append(line)
for i in loc_list:
    arcpy.Select_analysis("locality.shp", "AOI.shp", " \"LOCALITY\" = '%s'" % (i))

where 'LOCALITY' is a field in 'locality.shp'.

The error returned is "AOI.shp" already exists. I conceptually understand why this isn't working, but is there a way I can make a selection in python for the list of localities in 'locality.txt' and create a new shapefile/FC based on that selection. I only want one shapefile with all the selected localities in this shapefile. There are about 600 localities, so creating 600 individual shapefiles and merging them into one at the end of the script is something I would like to avoid.

I feel like there is a simple solution to this and I am overlooking something easy.

I'm also aware of list concatenation, but only recently aware, so haven't yet used it in any of my code.


texts:
I need to export color orthophoto map to gray scale geotiff file using free tools (like qgis) any idea?

I tried qgis preview but this does not produce new geotiff file


texts: mask raster with shapefile using gdalwarp/gdal

texts:
I am struggling to use gdalwarp with cutline option.
The following command keeps failing and I do not understand why :

gdalwarp  -cutline shapefile.shp  file.tiff newfile.tiff

  Warning 1: Self-intersection at or near point 1890.1578688362522 >211.13062073074084

  ERROR 1: Cutline is not valid after transformation

The produced file is indeed not valid being filled with NA ...
I tried with several shapefiles and tiff files but I could not make it works.

I use debian and gdal version 2.1.0


texts: Making Map Document (*.mxd) Locked/Uneditable?

texts:
I am wondering if its possible to create a locked map document?

I have a document that multiple users will be accessing and as a way of preventing them from making any accidental errors, am wondering if I can create a map document that is locked as to prevent any modifications.


texts: Automate Naming of Stormwater Network (Geometric network)

texts:
We have a geometric network representing our stormwater network. For anyone unfamiliar with stormwater drainage networks, they are composed of pipes and pits (simplified for this question). Pits are required wherever two pipes join or a change of direction is required. So in the geometric network, pits are junctions and pipes are edges.

The feature classes within the geometric network are a Point Feature for the pits and a polyline Feature for the pipes.

My task is attribute every pipe and pit with unique name in line with our Naming Policy. The name of the pipe or pit is dependent on its location within the network. Simplified naming policy is as follows:

Pits are given a number starting with 1 at the outlet (furthest downstream)
Pipes are given the number of their upstream pit

I would like to automatic this naming process. I have little experience in arcgis and python, but i have a background in Java. Does anyone have any advice as to how to achieve this? Manually naming the pits/pipes will take 4-8 weeks. I would rather spend 8 weeks learning ArcObjects etc. I’m guessing i need to represent the network as a graph in java, from there the naming process would be easy.

How do I extract the geometric network relationships into a graph?


texts: Restricting GeoServer to OpenLayers on localhost?

texts:
I'm publishing WMS-Layer with GeoServer that should be accessible only by an OpenLayers client running on the same machine.
The OpenLayers map is not public and I want to restrict also the WMS delivered by GeoServer so that it should only be visible to users accessing the layer through our OpenLayers client.

Since I'm new to GeoServer I actually don't see another way to achieve this, but there might be other and even better ways to restrict the access to the WMS-layers from GeoServer?

What should I try?


texts:
I am trying to merge 2 files, but the instruction gives this error:

"'ascii' codec can't encode character u'\xe0' in position 213: ordinal not in range(128) See log for more details"

what can it be?


texts:
I have been using gdal.Warp() in Python 3.6 for some time now and have just noticed the workingType option (see here for gdal.WarpOptions()). From its definition:

  workingType --- working type (gdal.GDT_Byte, etc...)

I take this to be the data type that Python converts the underlying data to, in order to complete the transformation/reprojection/warp.

My question is: What is the default working type used if one is not specified in the gdal.WarpOptions()?

If my data is in float32, does gdal keep this as float 32, or convert to another data type by default? I'm warping some fairly large files and am cautious of using too much memory etc.


texts: Several fields and GEOMETRY GetFeatueInfo from GeoServer

texts:
I have a table in PostGIS 2.0 with several fields of type GEOMETRY.
In the geometry_columns they are published in GeoServer 2.2 as layers. This is done as here - GeoServer docs, but immediately after the GetFeatueInfo stopped working in these layers.
Does anyone know whether you can configure the server to serve multiple fields of geometry.


texts: FME Use Output of reader in another reader

texts:
I have a ESRI geodatabase reader to read input table and i am using where clause to limit data based on date. write now i am using the user parameter to ask user latest date for input table. and then i am writing it to output table.

I want to check the maximum date in output table and based on that i want to give that maximum date as an input to input table reader so that every time my FME run user do not have to manually enter maximum date.

I have searched for this and now i have the maximum date as i have added the output layer as a reader and i am using Statistics Calculator to get maximum value.

Now How can i use that maximum value as an input to where clause of input table reader ? There seem to be no way of connecting anything to reader.


texts: MapInfo, SQL data from multiple tab files which have the same columns

texts:
New to the GIS forum so apologies in advance for any potential etiquette failures...

I'm working in MapInfo. I have multiple tables with the same structure. One of my columns is pipe diameter. I would like to do an SQL select for a particular pipe diameter or range of diameters i.e. 5m<= x <=10m across all the tables. Is this possible?

I am able to do it for one table and save the data into a new table file, then redo this using a Mapbasic script however, I then end up with lots of individual tables. I could do this and then copy and paste between the tabs however the data is very large so this is time consuming. Using the SQL tool, is it possible to select data from a particular column which reoccurs in multiple tables, and create a group selection across multiple tables which can be saved into a single table?

Feel free to ask any questions


texts: Better-matching bounding box by using collection

texts:
I am looking for a way to approximate a geographic region with a few bounding boxes. I know I can get a single box to cover the whole region, but some of my shapes are quite strange (gerrymandered, in some cases).

I'd like to know if there's a way to approximate a shape as a (small) collection of axis-aligned rectangles to make a tighter collection of bounding boxes.

For example, I know I can draw a rectangular box around a roughly "C" shaped region, but I can make a much better approximation with 3 rectangles like so:

    +--------------+
    |              |
    +----+---------+
    |    |
    |    |
    | +~~+--------+
    +-+~~+        |
      +-----------+

Constraints:

The rectangles should contain the entire original shape.
It's fine if the rectangles overlap.

Trade-off:

I'd really like to stay under 10 rectangles. (Maybe an increased penalty for each additional rectangle?)
I'd also like to have under 25% of the original shape's area as extraneously included by the rectangles, but if it takes too many shapes to do so, I'll go with fewer shapes.

Is there a recognized method for doing this to arbitrary shapes? To certain classes of shapes?


texts:
I am working on a project with two different teams, each working with a different GIS software package: QGIS and ArcGIS. I am trying to convert my QGIS layout to ArcGIS, so the team with ArcGIS can see exactly the same symbology as I do on my computer.

I am currently working with .qml files to save my layer styling, and I would like the other team to be able to visualize the same symbology for the layers we're both using. So I want to convert all my layer symbologies from .qml to .lyr.

Is it possible for them to use the .qml format in ArcGIS?

Is there a way to make all projects have the same symbols, without having to exchange RGB values, point style, line style, etc, etc...?


texts:
I am trying to convert Esri file gdb to PostGIS

I created a new database : 'newkelibia' with pgAdmin 4

using OSGeo4W shell I executed this :
 ogr2ogr, -f "PostgreSQL" PG:"host=localhost port=5432 user='postgres' password='1234' dbname='newkelibia'" C:\Users\ons\Documents\CityEngine\Default Workspace\new_project\models\KlibiaBD.gdb" -overwrite -progress --config PG_USE_COPY YES

and I got this :
ERROR 1: PostgreSQL driver doesn't currently support database creation.
Please create database with the `createdb' command.
ERROR 1: PostgreSQL driver failed to create

Did I mess any steps ?


texts:
I want to set up TileStache, but I'm new enough to the world of GIS to find the TileStache docs a bit impenetrable.

I've just about figured out that I need to run

tilestache-server.py -c tilestache.cfg

to start TileStache, and navigate to http://127.0.0.1:8080/ to see the results.

However, I'm unclear on what to put in tilestache.cfg.

Eventually I want to use GeoJSON, generated from a PostGIS database, but for now I'd just like to start with a simple example to make sure that TileStache is working.

There is a 'minimal sample configuration' in the docs, but it seems to rely on Mapnik - and refers to a style.xml file - do I need to install and configure Mapnik before I can use TileStache?

* UPDATE *

Okay, I've tried using a 'proxy' config as follows in tilestache.cfg:

{
  "cache":  {"name": "Test"},
  "layers":
  {
    "roads":
    {
      "provider":
      {
        "name": "proxy",
        "url": "http://tile.openstreetmap.org/{Z}/{X}/{Y}.png"
      }
    }
  }
}

And then restarting Tilestache and going to http://127.0.0.1:8080/, but I still see:

Bad path: "/". I was expecting something more like "/example/0/0/0.png"

If I go to http://127.0.0.1:8080/example/0/0/0.png then I get a 404 - any ideas?

* UPDATE 2 *

For the benefit of others: roads should be part of the URL, so visit http://127.0.0.1:8080/roads/0/0/0.png or http://127.0.0.1:8080/roads/preview.html to see tiles!


texts: Converting polygon in grid points for each grid of 100m*100m in PyQGIS?

texts:
I want to convert my polygon layers into point layer where I will have a point for each grid of 100m *100m grid for entire polygon layers using PyQGIS ?

I don't know any library function for grid creation and way to extract single point from the gird Point may be centroid of grid.

That's why I haven't provided any code snippet.


texts: Measuring line length in polygons with virtual layer qgis 3

texts:
What seemed to be an easy task turned out to be not that easy.
I want to measure line length in separate polygons.
My steps so far:

polygons to lines (because intersecting lines with polygons ends up in lineends out of the polygons get deleted, same behaviour like clip btw. not logic for me)
split lines with lines (polygon lines with actual line layer)
save intersection (because you cant apply a virtual layer on a temporary layer)
calculate line lengths (lines splitted need new line lengths)
Run the virtual layer with st_within
OR
select by location with "within", select splitted line layer with polygons one by one.

And here strange things happen: The lines, cutted with the lines from the polygon, still dont get correctly selected at the boarders. I get one half, but not both of the cutted lines.
Like, one part of the splitted line isnt within the polygon. But i cutted them with the boarders of the polygon and didnt move them.


texts: filterBounds() error in Google Earth Engine

texts:
I wanted to clip the ALOS elevation model to an ROI. This is the code I wrote.

var dem = ee.Image("JAXA/ALOS/AW3D30_V1_1");
var geometry =
    /* color: #ffc82d */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[1.3787373046875473, 49.38636777007675],
          [1.3787373046875473, 48.32016208111379],
          [3.5265644531250473, 48.32016208111379],
          [3.5265644531250473, 49.38636777007675]]], null, false);

var clipdem = dem.filterBounds(geometry)

But it gave me this error:

  dem.filterBounds is not a function

I tried casting it to an imageCollection before applying the function, but it doesn't clip the image - giving me the image for the entire world. What is causing it and is there a way to work around it?


texts: I am using the gIntersects function from the rgeos package to check whether a river layer (polyline) intersects with a species maps (polygons). However, it is taking ages to do the calculations, despite I have simplified the layers with gSimplify and added prepared=TRUE to the function formula. Is there any way to improve the speed of the process?


texts: Merge large numbers of ascii files .asc in ArcGIS

texts: I need to merge ~8k of .asc files into one big one. To convert each one into a raster first would take forever .. is there a possibility to convert multiple ascii files to raster in ArcGIS?


texts: ArcGIS Server: Execute Export Web Map Task with custom image service rendering rule

texts:
I am trying to use the out-of-the-box Print ==> Export Web Map Task for ArcGIS Server (10.1+) a web map containing an image service layer with a custom rendering rule. The web map json spec mentions the ImageService Layer has :

ArcGIS ImageService Layer

  {
    "id": "<webmapOperationalLayer_imageService>",
    "url": "<url1>,
    "token": "<tokenString1>,
    "title": "<title1>",
    "opacity": <opacity1>,
    "visibility": <true | false>,
    "minScale": <minScale1>,
    "maxScale": <maxScale1>,
    "noData": <no data value>,
    "format": "<jpgpng | png | png8 | png24 | jpg | bmp | gif | tiff>",
    "interpolation": "<RSP_BilinearInterpolation | RSP_CubicConvolution |
                            RSP_Majority | RSP_NearestNeighbor>",
    "compressionQuality": <compression quality>,
    "bandIds": [
      <bandId1>,
      <bandId2>
    ],
    "mosaicRule": "<mosaic rule 1>",
    "**renderingRule**": **-rendering rule 1-**
},

When the renderingRule is included, it receive a 404 error. When I do not set renderingRule, the task works fine.  I use the same rendering rule as when I call exportImage on the image service itself.

Has anybody seen an example of Export Web Map Task handling a custom renderingRule for an Image Service Layer?

Thanks in advance for your help.


texts: Setting variable from input map in GRASS python script?

texts:
I'm currently working on a GRASS script that takes maps from user input and computes them into a mapcalc expression. The issue i'm stucked with that when i ask for the map i can't set a variable name on it to use it on the mapcalc command. The code looks like this:

Asking for the map:

  #%option G_OPT_R_INPUT

  #% key: c

  #% key_desc: name

  #% gisprompt: old,raster,raster

  #% required: yes

  #% description: Effective cohesion (kPa)

  #%end

and trying to set a name on it:

  c = options['c']

I have previously imported parser from grass.script.


texts: arcpy.SubsetFeatures_ga output not saving in correct feature dataset

texts:
I am using SubsetFeatures_ga to create a random dataset within a file geodatabase and am trying to save the output feature class in a feature dataset but it keeps getting saved in the root of the file geodatabase. I have ArcMap 10.3.1 installed and am using Python 2.7 to run the code. I have already read this post which simulates the same problem but it was apparently solved by ESRI at ArcMap 10.1.

Here is the code:

import sys, string, os, arcpy
from arcpy import env
env.workspace = "J:/data/database/test.gdb/polygons/"
arcpy.CheckOutExtension("GeoStats")
in_features = "polys"
out_training_feature_class = "polys_train"
out_test_feature_class = ""
size_of_training_dataset = "10"
subset_size_units = "PERCENTAGE_OF_INPUT"
arcpy.SubsetFeatures_ga(in_features,out_training_feature_class,out_test_feature_class,size_of_training_dataset,subset_size_units)

Here is an image of the database structure:

And this is what it looks like after I run my script (you can see how 'polys_train' is output into the root of test.gdb)


texts: Saving QGIS form design basic functionality in .ui to modify it with QTDesigner

texts:
It is possible to save a classic Qgis form (a drag and drop form) in order to modify it with QtDesigner afterwards ?

I have created a form and I just want to change a list value into a toggle button but I don't find the way to save the form in .ui format.

I am using QGIS 3.4.5.


texts:
I am looking to create a 3D WebPortal for my client. My organisation has ArcGIS for Server so we would be looking for the easiest way to deploy using this (although open to other suggestions).

There is potentially a large amount of data...DTM/DSM, 3D buildings, Aerial Imagery and many vector datasets.

Have any of you deployed something like this before. I guess the key problem is the '3D part. Currently we supply a standalone appplication that uses ArcGIS Explorer, but they now want this across the organisation as a webportal. Or maybe served from our server to ArcGIS Explorer, if that can be done?

Any suggestions are much appreciated!

Regards


texts: What is this 'Run layer action' and how to use it?

texts:
Following on from this question: Run QGIS layer Action on all selected features, a possible solution was to create an action in Python which iterated through the layer's selected features.
If you right-click in the attribute table, there is a Run layer action which appears to be greyed out (Print information is an action):

What is this Run layer action option and how can it be used?


texts: QGIS 2.2.0 install via OSGeo4W overwrote/lost qGIS ver. 2.0

texts:
Just installed ver 2.2.0 Valmiera through the OSGeo4W advanced installer for a test drive.  It overwrote existing ver 2.0 in the process. Still want to have the previous to work from.

How can i get QGIS ver 2.0 back without problems in running the new release?


texts:
The top two layers in my table of contents are the most important. Streams and some point shapefiles. They overlay a raster. In layout mode I switched from portrait to landscape. Immediately the streams thinned so much they are barely visible and the points are so small they are invisible. (I know they are there because they blink when I turn the layer off and on). Correct locations though.

There were two similar topics online. One suggested checking the size and position in the properties to make sure it isn't negative. Done. The other was to delete the border. Done.

Unfortunately neither worked.

I'm on ArcMap 10.1 if that helps.


texts: Is it possible to edit related tables in ArcGIS Online?

texts:
I would like to make an editable web map where the user can create a feature then fill in information in a related table all in ArcGIS Online.  Is there any way to do  this?  So far, my research says no, but my organization is sold on the idea and would like to find a way to do this if there is one.

I am aware that it is possible to edit a related table using the ArcGIS Collector App as explained here, but we would like to work with ArcGIS Online on a desktop if possible.

Even if this is not currently possible, I'd like to know definitively so I can start exploring other solutions.


texts: Windows Server 2019 + Apache 2.4 + GeoServer 2.15.1 (HTTPS)

texts:
I have a server with Apache 2.4 on Windows server 2019. I´ve configured Apache to run as HTTPS with corresponding certificate.

I've installed GeoServer 2.15.1 for Windows that use Jetty how server.

I have an error of type

  Mixed Content: The page at 'https://localhost/' was loaded over HTTPS, but requested an insecure frame 'http://localhost/geoserver/UNO/wms?SERVICE=WMS&REQUEST=GetFeatureInfo&FORMAT=image/png&TRANSPARENT=TRUE&STYLES=&VERSION=1.3.0&LAYERS=secciones_mayor64&WIDTH=981&HEIGHT=799&CRS=EPSG:3857&BBOX=-919776.1230229893,5385864.567878975,-901029.941836075,5401132.864278918&INFO_FORMAT=text/html&QUERY_LAYERS=secciones_mayor64&FEATURE_COUNT=25&I=459&J=496'

  This request has been blocked; the content must be served over HTTPS.

I know that problem is because GeoServer runs on HTTP. How can I configure GeoServer or Jetty to use HTTPS?


texts: Generating shapefile of polylines using pyshp from a csv file and opening file in qgis

texts:
I tried the solution posted in:

Creating Polyline from GML using PyShp?,
Using PyShp to convert polygons in *.csv to *.shp files? and
Using PyShp to create polygon shapefiles?

to generate polylines from a csv file.

When I load the shapefile into QGIS, it is not displaying the lines.

Is this a known issue with QGIS?

I am using QGIS 2.8.2 on a Mac.

I can see that the coordinates are in the shapefile when I do this:

import shapefile
test = shapefile.Reader(shpFilePath)
shapes = test.shapes()

for i in xrange(0,len(shapes)):
    x = shapes[i]
    a = x.__geo_interface__
    print a['coordinates']

Is there another software we could use to view the shapefiles?


texts:
I need to retrieve some data from CSW server running GeoNetwork. I'm using this query but it doesn't work:

http://serverurl/geonetwork/srv/en/csw?
REQUEST=GetRecords&
SERVICE=CSW&
VERSION=2.0.2&
OUTPUTSCHEMA=http://www.isotc211.org/2005/gmd&
constraint=AnyText+EqualTo+Africa&
constraintLanguage=CQL_TEXT&
constraint_language_version=1.1.0&
typeNames=csw:Record&
RESULTTYPE=hits

Always returns to me the whole metadata (instead of 0; that's what I expect because there aren't any resources with name 'Africa').


texts: Ecw not supported with qgis-mapserv-fcgi but supported in QGIS Desktop

texts:
I have problem with qgis-mapserver version 2 and ECW.

Qgis Desktop works well with ECW (gdal support for ecw with the library ERDAS_ECWJP2_SDK-5.0) but not qgis-mapserver (not supported file in Apache error logs).

qgis-server must use the same gdal librarie that qgis desktop (QGIS desktop and server are all 2 installed on my computer (ubuntu 12.10). This is what I do not understand.

Note:

  I just watched the license to use the ERDAS_ECWJP2_SDK-5.0, including
  this paragraph: 1.9. "Server Application" means clustering software
  products That You Develop or maintain maintenance and are able of
  running on a Server To Provide an image decoding or processing
  capability. This License Agreement Does not Authorize You to use the
  SDK to create Server Applications.

  With this SDK for a "free" installation, or it is installed with a
  "Desktop Read-Only" license or with a license "Server Read-Only End
  User" but not the 2.

  Next time, I would read longer uses licenses, something I rarely do.


texts:
First time I used the leaflet. I want to create a popup on marker. On mouse-over of marker I want to get a popup message . my code is :

for (var i = 0; i <= result.length - 1; i++) {
    var px = result[i]['p1_x'];
    var py=result[i]['p1_y'];
    var map = this.state.maps[0];

   // create custom icon
    var firefoxIcon = L.icon({
         iconUrl:  LOCAL_BASE_URL + '/assets/heatmap/' +'/HARBOURFRONT/' + "/" + BASEMAPS + "/" +  'icon_sensor.png',
         iconSize: [15, 15], // size of the icon
   });

   console.log(" drawAllBaseMaps1 add icon inside map");
   // create marker object, pass custom icon as option, ad/d to map
   var marker = L.marker([px, py], {icon: firefoxIcon}).addTo(map);

   marker.on('mouseover', function(e) {
           //open popup;
          var popup = L.popup()
            .setLatLng(e.xy) //(assuming e.latlng returns the coordinates of the event)
            .setContent('
Hello world!
This is a nice popup.
')
            .openOn(map);
  });

It is not working.


texts: In QGIS Composer, is it possible to draw coordinates next to grid cross markers?

texts: I know how to add coordinates to the border of the map, both inside and outside. But I was wondering if there's a way to add coordinates next to grid markers, such that when I choose the option to have the cross markers, each cross has coordinates next to it.


texts: Is there a prj-file-defintion available for the Hölzel-projection? (http://dl.maptools.org/dl/libproj4/manual.pdf#page=58, page 58)?


texts: How to get the Area of two intersecting polygons on PostGIS?

texts:
I am trying to extract information from a base map based on overlaying polygons as shown in the picture.

The base map represents different groups (Names 1-8).

From the overlaying polygons (gray,names 1-7) I would like to know the following:

What is the area or proportion of the overlay (e.g.Polygon3= 40% (purple) in group 1, 60% in group2 (orange))

I tried to get there using:

SELECT
t1.name,
t2.group_name
FROM
public.overlay AS t1,
public.basemap AS t2
WHERE st_intersects(t1.geom,t2.geom)=TRUE
ORDER BY t1.name;

This gives me the underlying groups.

However, as some geometries of the base map are in the same group (same color) it would be good to have this summed up (e.g. polygon 4 intersects three geometries of the base map which are all belonging to the purple group).

I did this using R before but would prefer to do it in the database directly.

My preferred output would look something like this :

name    |group_name| coverage
1       |1         | 1
2       |2         | 1
3       |1         | 0.4
3       |2         | 0.6
4       |1         | 1
5       |1         | 0.65
5       |2         | 0.35
...


texts: Python Max Function in Field Calculator for Date values

texts:
I have 4 fields with Date Time values, and created a new field (Date type as well) which takes the most recent date from these 4 fields. I have tried using the Max python in the field calculator to populate this new field, but it does not use the most recent date from these 4 fields. Does this function not work with Datetime values, or is there another function for this that I am not able to find on these forums.

max([!field1!, !field2!, !field3!, !field4!])

I am on 10.5.  When I have two fields that have a date within the same month, it will not read the rest of the date correctly. (See picture attached)

When I run the following updated code it fails and says it is a syntax error.

With the current code, it reads the first number in the day to select the most recent date. (As shown below). It does not use the correct chronological format, even though the strptime format worked previously.


texts: Transform coordinates EPSG:32628 to EPSG:3857 Openlayers 3

texts:
I need to transform the example coordinates to a different projection in Openlayers3, the EPSG:32628 to EPSG:3857. I made a function to do it but the result is exactly the same numbers.

function plsHelp(){
    var pt = [459123.1209, 3108334.7701];
    var scr = new ol.proj.Projection('EPSG:32628');
    var dest = new ol.proj.Projection('EPSG:3857');
    var test= ol.proj.transform(pt,scr,dest);
    console.log(test);
}

What is wrong with the function that is preventing the tranformation of the coordinates from EPSG:32628 to EPSG:3857 in Openlayers 3?


texts:
I am getting started with ArcGIS JavaScript 4. I want to create a layout with header menu, left pane and right pane with map. All the ESRI examples seem to refer to JavaScript 3. In the code below when I try to use this with JavaScript 4 API I get a 404.

Is there a difference in the way layouts are created in JavaScript 4 API?

Is it still with Dojo?

In the head tag I have :

  <head>
<link rel="stylesheet" href="https://js.arcgis.com/4.13/esri/themes/light/main.css">
 <script>var dojoConfig = { parseOnLoad: true };</script>
<script src="https://js.arcgis.com/4.13/"></script>

then in script tag after head is closed:

   <script>
  dojo.require("dijit.layout.BorderContainer");
  dojo.require("dijit.layout.ContentPane");
  dojo.require("dijit.layout.TabContainer");
  dojo.require("esri.map");
  dojo.require("esri.arcgis.utils");
  dojo.require("esri.IdentityManager");

require([
  "esri/Map",
  "esri/views/MapView"
], function(Map, MapView) {

var map = new Map({
  basemap: "topo-vector"
});

var view = new MapView({
  container: "viewDiv",
  map: map,
  center: [174.80500, -36.82700], // longitude, latitude
  zoom: 11
 });
});

This results in This code results in a 404 error for https://js.arcgis.com/4.13/dojo/dijit.layout.BorderContainer.js net::ERR_ABORTED 404


texts: Calculating length of polyline using ArcPy?

texts:
I have a polyline shapefile that has only one polyline.

How do I calculate the length of that polyline using ArcPy?


texts: Performing relate in ModelBuilder for one to many relationship?

texts: I'm not sure how to do this as I don't' think there are any ModelBuilder tools for this task.  I'm not real savvy with Python but basically I want to relate field1 in table A to field2 in table B with a one to many relationship.  So I want all records in table B that match Table A.


texts: How to reference a specific layer in a dataframe with Arcpy

texts:
I am trying to reference a Layer named "Tim" in the dataframe "Layers" in my .mxd. It is not working. I have looked at the ArcGIS online help and all else in this forum. Can't get rid of this Runtime Error!

>>> mxd = arcpy.mapping.MapDocument("CURRENT")
>>> df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
>>> layer = arcpy.mapping.ListLayers(mxd,"Tim",df)[0]
Runtime error <type 'exceptions.IndexError'>: list index out of range

Any insight would be greatly appreciated!!

So this is my entire script:

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "layers")[0]
targetGroupLayer = arcpy.mapping.ListLayers(mxd, "Tim", df)[0]
addLayer = arcpy.mapping.Layer(r"MyLayer.lyr")
arcpy.mapping.AddLayerToGroup(df, targetGroupLayer, addLayer, "BOTTOM")

this also give me the error.

When i change

targetGroupLayer = arcpy.mapping.ListLayers(mxd, "Tim", df)[0]

to this

targetGroupLayer = arcpy.mapping.ListLayers(mxd, "", df)[0]

it works but i need it do reference a specific group layer because there will be many of them and eventually i will pass in a variable for the targetGroupLayer so i can bring in certain layers into that group.

My most recent as Roy suggested.

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
for lyr in arcpy.mapping.ListLayers(mxd, "",df):
    if lyr.name == "Tim":
        addLayer = arcpy.mapping.Layer(r"C:\Users\T\Desktop\MyLayer.lyr")
        arcpy.mapping.AddLayerToGroup(df, lyr, addLayer, "BOTTOM")

No luck still!

Thanks for all the suggestion also!! I will keep trying!

** This script is working now!!* Thanks Roy!


texts:
I have a folder with .adf raster data that I'm trying to import into ArcMap 10.6, but I can't figure out how. I'm trying to do menu File → Add Data → Add Data, but the files are not displayed. Nor do they show up under "Contents" in ArcCatalog, but in preview the raster is displayed.

Online instructions talk about "import", but I don't see any menu item that's called this.

How can I import them?

Data is from this website: 50 m Arc GIS Grid
http://www.soest.hawaii.edu/hmrg/multibeam/bathymetry.php


texts:
I have used symbol levels in my current project to define the hierarchy of the different roads. I have tried to make the symbol levels the same in order to achieve the desired look, with no luck. The roads are a minor element in the map that I am making and will all be symbolized the same except for width.

Shown in this screenshot:

(Colors have been added for clarification. The green are US Routes, orange are State Routes and white are local roads.)

If they are all white, is there a way in QGIS to have them look joined at their intersections, similar to all intersections of the local roads?

This is a tough one to explain, so if clarification is needed, please let me know.


texts: Accessing non-PostGIS Postgres instances from QGIS

texts:
Does anyone know if it's possible to query tables that are regular Postgres and not PostGIS with QGIS? I have a PostGIS database but want to access tables that aren't spacial - much like reading in a text delimited file.

I've poked around the plugins but to no avail


texts: How do you draw an adjacent polygon in JOSM?

texts:
In JOSM, how do you draw an polygon that is adjacent to another polygon without having to retrace the shared boundary? I'm looking for the same functionality that is found in Auto-Complete Polygon tool found in Esri software: http://desktop.arcgis.com/en/arcmap/10.3/manage-data/creating-new-features/creating-adjoining-polygons.htm

What I've tried to do is draw the new polygon boundary (that doesn't have a shared boundary). Then use the Split Adjacent Ways tool. If I select the new boundary, the split boundary, and then combine them, it works but I break the other polygon. I also tried to Copy/Paste or Duplicating the split boundary, so there are two ways. I then would select the new Way to combine. But that doesn't work because the copy/paste/duplicate doesn't put the geometry in the same place.


texts: PgRouting routing algorithm from start to end via nodes with restrictions

texts:
I have the following situation:

I want to find the shortest route from A to B through nodes (C,D,E, etc.) with turn directions and road type filtering.

Is there a way to accomplish this using pgRouting? Has anyone done this? I haven't begun developing anything yet, I'm just trying to find out if pgRouting is the way to go.

Thanks.

Update: The algorithm should determine the order of nodes (C,D,E, etc.)


texts: Wrong results joining places to create flow map

texts:
I have to do a flow map (home->working place) for my area of interest (28
french cities).

I get a database .dbf created by a french institute. I imported it into
postgres/postgis.
One row = 1 people.

My table is composed by several fields :

id (primary key)
codgeo_res: residence (home) city 's id
x_resi: residence (home) city 's longitude(centroid)
y_resi : residence (home) city 's latitude(centroid)
codgeo_tr : work place city 's id
x_tr: work city 's longitude(centroid)
y_tr :work city 's latitude(centroid)
ipondi : poids de l'individu (weighting value assigned to each individual du to the sampling)
geom : residence (home) city 's ponctual geometry (centroid)

So my goal is to create a new table in order 1 row = 1 flow (city A toward
city B)

I have to generate useful  geometries (points + lines arrow). Moreover, have
to apply a SUM on the ipondi field and a GROUP BY on "codgeo_res" field.

Could you help me in order to translate this
objective to functional sql queries?

Up to now, i tried to run this query throw qgis interface called "virtual
layer creation" :

------------------------------------------------------
Code:

select
st_makeline(st_makepoint(entrants_only.x_resi,entrants_only.y_resi),st_makepoint(entrants_only.x_tr,entrants_only.y_tr))as
geometry,sum(entrants_only.ipondi) as sumipondi, entrants_only.codgeo_res as
codgeo_res, entrants.codgeo_tr as codgeo_tr

from entrants_only
group by entrants_only.codgeo_res
--------------------------------------------------------

But i get this error message :

---------------------------------------------------------
query preparation error on PRAGMA table_info(_tview): no such column
entrants_only.geom
---------------------------------------------------------

Here's a screenshot representing the error message.

The problem was due to the datatype of my XY fields (string). So I did a cast STRING to BIGINT (those fields are storing integer EPSG 2154 Lambert93).
So now, the query works. I attach a screenshot.

Indeed, the postgis temp linear layer is created.
But I have always a problem regarding the GROUP BY. I want to group by on one field only (codgeo_res), but I can't generate the query without put all
fields present in SELECT into GROUP BY!

This is my current query :

create table temp as select
st_makeline(st_makepoint(entrants_v2.x_resi,entrants_v2.y_resi),st_makepoint(entrants_v2.x_tr,entrants_v2.y_tr))as
geometry,sum(entrants_v2.ipondi) as sumipondi, entrants_v2.codgeo_res as
codgeo_res, entrants_v2.codgeo_tr as codgeo_tr

from entrants_v2
group by
entrants_v2.codgeo_res,entrants_v2.x_resi,entrants_v2.y_resi,entrants_v2.x_tr,entrants_v2.y_tr,entrants_v2.ipondi,entrants_v2.codgeo_tr

Is there a way to assign one field into group by with PostGIS?

My result is false.


texts:
I'm trying to develop a WFS-T on a web page using the Boundless WebAppBuilder SDK with QGIS.  I have an HTML page referencing a WFS coming from GeoServer with an MSSQL source.  When I try to draw a new feature, I get an error from Firebug saying:

"No such feature type http://(my ip):(port)/geoserver/DRTEST:polygontest"

This is the section of my JavaScript with the parameters:

var lyr_polygontest = new ol.layer.Vector({
     opacity: 1.0,
     source: wfsSource_polygontest,
     style: style_polygontest,
     selectedStyle: selectionStyle_polygontest,
     title: "polygontest",
     id: "polygontest20161121124618872",
     wfsInfo: {
         featureNS: 'http://(ip:port)/geoserver/DRTEST', //http://(ip:port)/geoserver/DRTEST
         typeName: 'DRTEST:polygontest',
         featureType: 'polygontest',
         geometryType: 'Polygon',
         geometryName: 'the_geom',
         url: 'http://(my ip):(port)/geoserver/DRTEST/ows?service=WFS&version=1.0.0&request=GetFeature&outputFormat=text/javascript'
    },
    isWFST:true,
    filters: [],
    timeInfo: null,
    isSelectable: true,
    popupInfo: "

Attribute	Value
id	[id]

"
});

Any suggestions on what parameters need to be changed?  I'm still trying to figure out if the problem is if GeoServer is only looking for the namespace instead of an actual HTTP request.


texts: How to get right reflectance from Sentinel-2 LIC images?

texts:
I found that the range of DN from Sentinel-2 images is 0 ~ more than 10000. How can I get reflectance(0~1)?

This picture showed the way to get reflectance. But in this way, reflectance will over 1. Reflectance over 1 means what? Abnormal value?


texts: tiff to shapefile using GDAL - slow speed

texts:
I am converting a 500MB .tif file to .shp file using GDAL.

I am using the following:

gdal_polygonize.py .tiff file -f “ESRI Shapefile” out.shp

However, I abandoned the first attempt as in 12hrs it has 2% progress.

I split the .tif file into 20 pieces and the first chunk took 5hrs.

Is there a more efficient means of doing this conversion?

The ultimate aim is to convert .tif to GeoJSON. I am very familiar with the .shp to .json conversion and I am thus using .tif to .shp as a stepping stone.

The file in question is vegetation data from http://rmgsc.cr.usgs.gov/outgoing/ecosystems/Africadata/

The file in question - af_labeled_ecosys.zip

The intention is display the vegetation data.


texts: I have a list of coordinates with which i can create a path and I want to export this list as shp. I have "shapefiles" package but I didn't find the way to do this. Is there a way to do this with this package or I need something else?I use  R


texts: Unable to merge CSV and CSVT with shapefile using QGIS due to "text" fields?

texts:
I know the question is all over this site, unfortunately I cannot find a solution. I have this CSV file and this CSVT file that I want to merge with this shapefile. I've tried to format the .csv using all software mentioned in the forum including Excel, Numbers, Open Office, Libre Office, Text Edit and Sublime text.

The .csv file is imported almost correctly as you can see from the pictures, but the fields are wrong. Therefore the merge function does not work, because the csv fields are still seen as "text".

Can someone give me a hand on this?


texts: Placing multiple data frames in a single map layout in ArcGIS Pro

texts:
I am trying to prepare a layout for my study area by method showing country to state to city in ArcGIS Pro, but all the data frames are showing the same area and when I try to change the data, the data in all the data frames is getting changed.
How do I add different data frames in a single layout in ArcGIS Pro.


texts:
I am trying to download global GeoTIFFs following the instructions here: https://git.wur.nl/isric/soilgrids/soilgrids.notebooks/-/blob/master/markdown/webdav_from_Python.md
However, I get the error "Received a NULL pointer". Here is the code I used:
import os
from osgeo import gdal,ogr,osr

outdir = os.path.join(r"C:\Users\myusername", "SoilGrids")
outpath = os.path.join(outdir, "sgtest.tif")

kwargs = {'format': 'GTiff', 'creationOptions': ["TILED=YES", "COMPRESS=DEFLATE", "PREDICTOR=2", "BIGTIFF=YES"]}

ds = gdal.Translate(outpath,
                    "/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/ocs/ocs_0-30cm_mean.vrt",
                    **kwargs)
del ds

I was able to use gdal translate to download a sample Landsat file, so my gdal installation (inside Anaconda) works with vsicurl in some circumstances. I'm not sure what I'm doing wrong when it comes to SoilGrids. This was what worked for my Landsat test:
ds = gdal.Translate(outpath, "/vsicurl/http://landsat-pds.s3.amazonaws.com/L8/001/003/LC80010032014272LGN00/LC80010032014272LGN00_B1.TIF")


texts: Removing small polygon gaps in Shapely polygon?

texts:
I know this issue has been addressed before with different software, but how do I do it with Shapely?

I have a polygon with a very small gap on the interior of the shape:

Here is a close up:

The polygon is valid so I can't use the buffer trick, and simplifying the polygon also doesn't work.


texts: Is there a way to convert/recover a feature class that was edited in a higher ESRI license to work with it in lower license?

texts:
I have a geodatabase where one of the feature classes was edited using an Advanced license, but I work in a Basic license and when I go to edit the feature class I am unable to.

Are there any tricks to copying this feature class into a new one in order to be able to edit the data in Basic instead of Advanced?  I tried copy/pasting the feature class within the GDB but got an ESRI license error.  Short of building a new feature class from scratch with the exact same schema and then copy/pasting the records over, is there any other way to get around the license limitation to just edit the data's spatial?


texts: How to fix cql_filter problem with unicode characters in GeoServer 2.17.2?

texts:
I am using GeoServer 2.17.2 and OpenLayers v6.4.3, I can not get result when my cql_filter query is consist of unicode characters, I found out this issue can solve by add a code in tomcat server.xml, But I didn't install tomcat therefore I can not use that solution moreover I tried to use encodeURIComponent() but it didn't work also. my query in cql_filter is "province = زنجان" for example, as well as this issue cause to I get cors() error:
<ServiceExceptionReport xmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.3.0" xsi:schemaLocation="http://www.opengis.net/ogc http://localhost:8080/geoserver/schemas/wms/1.3.0/exceptions_1_3_0.xsd">
<ServiceException> Could not parse CQL filter list. Lexical error at line 1, column 12. Encountered: "\u0632" (1586), after : "" Parsing : province = زنجان. </ServiceException>
</ServiceExceptionReport>

How I can fix this problem?
     var cqlWmsLayers=
     new ol.layer.Tile({
         title:`${layerName}`,
           visible: true,
      source: new ol.source.TileWMS({
       url: `http://localhost:8080/geoserver/${workspaceName}/wms`,
       params: {'LAYERS': `${workspaceName}:${layerName}`,'cql_filter': "province = زنجان" },
 ratio: 1,
        serverType: 'geoserver',

OS: Debian 10 x64


texts:
I joined some fields to my tables and now I need a tool or something to change that field's name.

How can I do this in ModelBuilder?

I didn't find any tool to do that.


texts: Building expression that involves fields in Layer Properties (size, rotation) using QGIS?

texts:
I need to build an expression in QGIS that involves fields like Size, Rotation or Offset.

Is there any method how I could involve these fields when I am building the expression?

So I have a Point layer and size of points is set to eg. 3 (Map unit) and I have this Column MERITKO in Attribute table that I need to change the size of each point so what I need is to build the expression: Size * MERITKO

How can I do that?


texts: Specifying coordinate reference system (CRS) for scalebar that is different to project CRS?

texts:
Is it possible to specify a coordinate reference system (CRS) for the scalebar that is different to the project CRS?

For example, my project CRS uses geographic coordinates (lat/long), the reason being that I need a lat/long graticule on the map. However I would like the scalebar to be expressed in kilometres (i.e. projected distance units, rather than angular measures).

An alternative solution would be for QGIS to have the capability to specify the CRS of the map grid / graticule, as in ArcMap.

This is not simply a question about how to select scalebar units (e.g. miles / kilometres / metres).


texts:
I'm using QGIS 2.18.20. On MacOS High Sierra 10.13.5
Only recently I haven't been able to find uninstalled plugins. The settings show the repository is "unavailable".
I'm not using a proxy. I've tried using http:// instead of https:// but it comes up with SSL error.

I've manually downloaded the desired plugins (Animove and Point Sampling Tool) but I can't find the .qgis folder on my mac (perhaps this is a lack of knowledge about macs so please advise on this mac users!) so I'm not sure how to go about manually adding the unzipped files


texts: Get a point on a polyline from end points given distance along poly line

texts:
I have start point and end point of a polyline. How can I get a point on that polyline from end points specified by the given distance.

How can I get this by using arcpy provided that this script should work on ArcView without importing any modules (Like this answer). I'm trying not to hard code values.

Related Question:

How to create a point along a line given distance
Interpolating points along line in Python for ArcGIS Desktop?


texts:
I'd like to quickly inspect the features in a feature class by double clicking the row in the table view and looking at the surrounding base layer.

This works, but the feature fills the screen, so I have to zoom out for each feature I inspect.

Is there a way to set the zoom layer, or another method for quickly inspecting the area around many features in a feature class?

Thanks!


texts: Joined attributes not available in Feature Service

texts:
I have a joined feature class (PROPERTY) and table (OG_LEASE). When I publish this as part of a map service and feature service, attributes are only available in the map service.

Map service fields:

PROPERTY.OBJECTID
PROPERTY.TAX_MAP_NUMBER
PROPERTY.OWNER
OG_LEASE.OBJECTID
OG_LEASE.COUNTY
OG_LEASE.TERM
(etc.)

Feature service fields:

OBJECTID

How can I ensure that the feature service has all the joined attributes (as well as the original feature class attributes!) published? The only warning that came up when publishing was that the join field didn't have an attribute index.

(I did look at Problems displaying joined attribute data through ArcGIS Server but it's not directly related to my question -- both my feature class and table are in the same ArcSDE database. I'm using ArcGIS 10.2.)


texts: Why is a Merge process via ArcObjects extremely slow?

texts:
In one of my arcobjects-based GIS tools I have a process which takes an array of shapefile paths and uses them as inputs with the ESRI.ArcGIS.DataManagement.Merge tool. I run it via a Geoprocessor instance. I use this setup for various other gp tasks (such as Get Count and Calculate Field), but I have found that a Merge process takes a much longer period to process than if I had run it in ArcMap or ArcCatalog. For example, a process which I have clocked at around 3 hours took over 27 to complete when performed with this method in ArcObjects.

Has anyone experienced this before, and are there any known methods for expediting the process?

EDIT:
To be more explicit about the workflow taking place here:

Create a backgroundworker object to do this in a separate thread to avoid UI lockup.
List shapefile paths using Directory.GetFiles and append them to a semicolon-delimited path string for the merge operation
Perform a Get Count on each Shapefile
Use the path string to perform a merge operation.

I've tried using IGeoProcessor and Geoprocessor to see if there is any noticeable increase in performance, but to no avail. There are no tight loops involved (unless you consider rapid-fire Get Counts), and I certainly have no use for an EditSession. For some reason, however, the process just plods along at an amazingly slow pace.

It is worth noting that I am working with unreasonably large data sets (the road network of North America, for example), so I doubt any method that makes use of RAM as a workspace will be appropriate.


texts: Keeping built-in colour palette with Google Earth Engine?

texts:
How can I make a built-in palette in EE apply to all bands of an Image(), not just a specific one?

Context: I am using Google EarthEngine to display data from the Cropland Data Layer (CDL), which is a cropmap. Typical CDL images come with a band cropland, which has a built-in colour palette (showing the same colour scheme as the original CDL).

The problem: Now when I do a conversion from "Image Collection with many images, each including cropland" to "Image with multiple cropland bands, where the name of the band represents the year", I loose the colour attributes! Ironically, if I were to rename a specific band to cropland when visualising, the palette will come back.

This seems to  suggest that the palette is associated to a specific band name (indeed, it shows up in properties as: cropland_class_names).

How can I make the palette apply to any band in the Image, with arbitrary band name?

Code:

var CDL_12 = ee.Image('USDA/NASS/CDL/2012')
var CDL_13 = ee.Image('USDA/NASS/CDL/2013')
var CDL_12_13 = CDL_12.select("cropland").rename("CDL_2012").
  addBands(CDL_13.select("cropland").rename("CDL_2013"))

print(CDL_12_13, "CDL_12_13")

First one (untouched) will show the right colours. Second (band renamed), won;t, while third (renamed band name dback to cropland) will show right again:

Map.addLayer(CDL_12.select("cropland"), {}, "CDL 2012 raw")
Map.addLayer(CDL_12_13.select("CDL_2012"), {}, "CDL 2012 from Image")
Map.addLayer(CDL_12_13.select("CDL_2012").rename("cropland"), {}, "CDL 2012 from Image: works")


texts: Geoserver SLD: Dynamic Variable Substitution

texts:
I'm facing a problem with Geoserver and its styling mechanism that I'm not sure whether it's possible to be solved or not. But first, a quick context...

Context

We want to use Geoserver and OpenLayers to display a map in our application. In this app, we will be able to choose a Brazilian state from a select field in order to see every city from this state. Also, consider that each city will be filled with a specific color.

The colors will be sent through a WMS Request instead of computed or defined inside the SLD, for example, using data from our layer dataset. That's the problem... One more thing to take into account is that some states in Brazil have more than 400 cities.

Solving the Problem?

The first approach we tried consists of creating a specific Rule for each city. However it turned out to be unfeasible because we have more than 300 cities to be displayed in some states. Imagine how big this SLD file would be!

The second approach we thought about was to dinamically get each color from the QueryString. For example, when displaying the city named "Salvador", it must search the QueryString for the parameter color-salvador. But we don't know how to it. Actually, it seems impossible. Am I right?

Take a look at the SLD below.

<Rule>
    <Name>City Rule</Name>
    <Title>City</Title>
    <PolygonSymbolizer>
        <Fill>
            <CssParameter name="fill">
                #<ogc:Function name="env">

                    <!-- Problem here. We need to define this value dinamically -->
                    <ogc:Literal>salvador</ogc:Literal>
                    <ogc:Literal>234545</ogc:Literal>
                </ogc:Function>
            </CssParameter>
        </Fill>
    </PolygonSymbolizer>
</Rule>

What we would like to do is instead of defining the literal salvador, inside the env function, we would like to use a Property like city-name that comes from our layer dataset. Something like this...

<Rule>
    <Name>City Rule</Name>
    <Title>City</Title>
    <PolygonSymbolizer>
        <Fill>
            <CssParameter name="fill">
                #<ogc:Function name="env">
                    <ogc:Literal>
                        <ogc:PropertyName>city_name</ogc:PropertyName>-color
                    </ogc:Literal>
                    <ogc:Literal>234545</ogc:Literal>
                </ogc:Function>
            </CssParameter>
        </Fill>
    </PolygonSymbolizer>
</Rule>


texts: API that can give the latitude and longitude ranges from the address

texts:
Is there any API available that can get the latitude and longitude boundaries from the address.

For example if we provide: 12380 Sunnyvale avenue, CA, 90005

Can we get the latitude minimum and maximum, longitude minimum and maximum for that address?


texts: Get Number of Rows, ArcObjects (VBA 6.0)

texts:
I am trying to count the number of rows in a feature class attribute table. I have declared a feature cursor which i want to use later in the code to do some processing and field value calculation. (In the code below pFClass was declared to the IFeatureClass interface)

Dim pCursor As IFeatureCursor
pCursor = pFClass.Search(Nothing, False)

Dim pFeature As IFeature
Set pFeature = pFCursor.NextFeature

Dim LngCounter As Long

Do Until pFeature Is Nothing
    LngCounter = LngCounter + 1
    Set pFeature = pCursor.NextFeature
LOOP

When i run this code to get the number of items in the table i use a msgbox to print the number to screen which is accurate. But then i Cannot reuse my pFeature because it is now nothing.

Other permutations i have tried include declaring 2 cursors and 2 IFeature objects and then using the first to get the total count and then the second to do subsequent calculations of the attribute table but to no avail. I looked at the ITable.rowCount but could not get that to work for me. I'd really appreciate some help on this.

Thanks


texts: How to query attributes cross-table in QGIS?

texts:
is it possible in QGIS to query a layer with data from another layer?
And if so ... how is the syntax for that?

Usecase: I have a lot of image-position points and a lot of monitoring sites.
I would like to set the monitoring area layer, so that it only shows areas where we also have images for this year.
I then use the extent of the monitoring area as the coverage layer to render my image-position maps with the Atlas plugin.
Without querying the areas to reduce them, i get a lot unneccesary maps with no image-positions, so i want to get rid of those areas which do not match (and of course not manually)

So, is there kind of a query to find out if a value in table1 column xy exists also in table2 column xy ?


texts: rt_raster_gdal_warp error while importing raster data into PostGIS

texts:
I am trying to import raster data into PostGIS. Documentation tells me to run raster2pgsql, which I do with the following command:

raster2pgsql -I -C -M -F -s 910001 test.tif mytable | psql -d postgis_test -U postgres

Here, SRID 910001 is self-defined to match the data that I have (check below). Now, when I run this, I get the following error:

FEHLER:  rt_raster_gdal_warp: Could not create GDAL transformation object for output dataset creation
KONTEXT:  SQL-Funktion Â»st_transformÂ« Anweisung 1

(FEHLER means ERROR and KONTEXT means CONTEXT). Googling the issue leads to http://126kr.com/article/5byaz3d0uvp, which suggests checking if PROJ is defined and PROJSO set correctly. So I have set PROJSO="libproj-9.dll" and SELECT postgis_full_version() yields

POSTGIS="2.3.1 r15264" GEOS="3.6.0-CAPI-1.10.0 r4265" PROJ="Rel. 4.9.1, 04 March 2015" GDAL="GDAL 2.1.1, released 2016/07/07" LIBXML="2.7.8" LIBJSON="0.12" RASTER

I suppose that is fine. I have defined the projection as follows with WKT string

PROJCS["WRF experiment 1",GEOGCS["WRF Sphere",DATUM["unknown",SPHEROID["WRF  Sphere",6370000,0]],PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",30],
PARAMETER["standard_parallel_2",60],
PARAMETER["latitude_of_origin",47.83],
PARAMETER["central_meridian",11.01],
PARAMETER["false_easting",36.5],
PARAMETER["false_northing",-30],
UNIT["unknown",15000]]

and PROJ string

+proj=lcc +lat_1=30 +lat_2=60 +lat_0=47.83 +lon_0=11.01 +x_0=547500 +y_0=-450000 +a=6370000 +b=6370000 +to_meter=15000

(I generated the WKT string from the PROJ one using gdalsrsinfo). Visual inspection with QGIS shows that the PROJ string works.
I am on Windows using Postgis 9.6. I have no idea on how to further debug this. Is there a way to get more verbose output? I can provide my data if it helps.


texts: Reference to layergroup object in Leaflet?

texts:
I want to build my custom layers control with JQuery.

My JS is like:

$(document).ready(function() {
    var map = L.map('map', {
                    center: [39.73, -104.99],
                    zoom: 10
    });

    var littleton = L.marker([39.61, -105.02]).bindPopup('This is Littleton, CO.'),
        denver    = L.marker([39.74, -104.99]).bindPopup('This is Denver, CO.'),
        aurora    = L.marker([39.73, -104.8]).bindPopup('This is Aurora, CO.'),
        golden    = L.marker([39.77, -105.23]).bindPopup('This is Golden, CO.');

    var cities = L.layerGroup([littleton, denver, aurora, golden]);
    cities.addTo(map);

    // use jQuery to listen for checkbox change event
    $('div#layercontrol input[type="checkbox"]').on('change', function() {
        var checkbox = $(this);
        var layer = checkbox.data("layer");

        // toggle the layer
        if (checkbox.prop('checked')) {
            map.addLayer(layer);
        } else {
            map.removeLayer(layer);
        }
    })
});

My html is like:

        <input type="checkbox" data-layer="cities">Cities

Here the map.addLayer and map.removeLayer failed due to Error: The provided object is not a Layer..

How can the layergroup be referenced by the "data-layer" attribute associated with the checkbox for it?


texts: ogr2ogr: 'ERROR 1: Wrong field type for ID' with primary key

texts:
I have a shapefile with the following structure:

I'm trying using ogr2ogr to import the above shapefile inside my PostGIS database with the following code:

ogr2ogr -append -f "PostgreSQL" PG:"host=localhost dbname=vector_data
user=postgres password=postgres"
grid.shp -a_srs EPSG:4326 -skipfailures -lco FID=ID
-lco COLUMN_TYPES=ID=VARCHAR -lco GEOMETRY_NAME=geom

The problem is when I run the above code, I get the following error message:

ERROR 1: Wrong field type for ID

I've understood the problem is because my ID field is a String and not an Integer, however I put COLUMN_TYPES=ID=VARCHAR on my script and ogr2ogr was not able to understand it.

I tried to import this shapefile using no field as FID and it worked, but I would like to import it in my database using ID as a primary key even it is a String.

Is there a way of doing it?

The shapefile is here: http://www.dpi.inpe.br/newcatalog/tmp/grades/grade_cbers4_mux.zip


texts: Round numbers in concatenated labels in QGIS

texts:
I've created some multi-line labels in QGIS 2.18.12 using the expression:

 "PIPE_DIA"  || 'mm' ||  ' ' || "PIPE_MATRL" ||  '\n'  ||  'US Depth: ' ||   ( "START_COVELEV" - "START_INVELEV")  || '\n' ||  'DS Depth: ' ||  ( "END_COVELEV" - "END_INVELEV" )

All ELEV fields are double with precision to 3 decimal places.

The output however randomly spits out an extra 11 decimal places like this:

I've gone into label properties and set decimal places to 2 but it does nothing. I've also tried using the concat( function instead but it returns the same error.

Anyone have any ideas or work arounds?


texts:
I have several sets of choropleth maps showing ranges of poverty by census tract, each using a different poverty measure. All of the maps use the same categorical ranges to define poverty status.

The maps clearly differ visually in the distribution of poverty, but I am having trouble finding a way to provide a statistical comparison between maps. I have tried using contingency tables and the comparison of chi squared values - however, this really only compares the counts of areas that fall into each category, but provides no indication of spatial coincidence.

Are there any statistical measures to compare these types of maps both numerically and spatially?

Edit: Examples of maps are below


texts:
I'm new to Python and I'm trying to create a script which adds fields to a feature class. The provided script is running well, but I want to create an additional parameter, which asks for the required field length of the fields to be created. The field length is only important if the field type "TEXT" is selected. The script should only take the additional parameter into account if the user selects "TEXT" under the field type. I think I need to put an if...else statement, but I don't know how do that.

    import arcpy
    # input data set (feature class, shapefile)
    inputFC = arcpy.GetParameterAsText(0)
    # list of field names
    inputString =arcpy.GetParameterAsText(1)
    fieldList =inputString.split(";")
    fieldType = arcpy.GetParameterAsText(2)

    for name in fieldList:
        arcpy.AddField_management(inputFC , name, fieldType)
        arcpy.AddMessage ("field created: " + name)

    arcpy.AddMessage ("Calculation complete")


texts:
How can I orthorectify an oblique image of a flat landscape when I know the conditions of the photo (coords, height, direction, lens, etc)?

I would like to perform this kind of stuff from GRASS orthorectification but in R.


texts: Is there a pdf form of documentation for Survey Connect?

texts: I found a paper for Updating field types in Survey123 and wonder if it comes from a larger volume .pdf source with more information documented.


texts: Defining Feature Set parameter for ArcGIS Python Toolbox tool?

texts:
I'm trying to create an ArcGIS Python Toolbox tool which uses a polygon for input.

I want the user to be able to create the polygon when the tool is run.

I know how to do this using both a ModelBuilder model A quick tour of using Feature Set and Record Set and for a Python Script tool in a normal toolbox. This is done using a Feature Set variable by setting the schema to a polygon feature class.

I know from here Defining parameter data types in a Python toolbox that the datatype for Feature Set in a python toolbox is GPFeatureRecordSetLayer.

If I define the parameter as below it is not possible to create a polygon at run time. I’m almost 100% sure that this is because I have not set a schema. But does anybody know how to set the schema or know if there is some other way to interactively create a polygon and use as input for a Python Toolbox tool?

InteractiveParameter = arcpy.Parameter(
displayName="Interactive polygon",
name="toolname",
datatype="GPFeatureRecordSetLayer",
parameterType="Required",
direction="Input",
enabled = True)


texts: In an attribute table, how should I go about adding a Unique ID into a field? For example, I have a Water Meter points layer with an attribute table but I need a field that list them numerically with WM attached to the numbers. I would like to have a field displaying WM1001, WM1002, etc. Is there any way to accomplish this? Thanks.


texts: Area calculation of mapped categories in ArcGIS Desktop?

texts:
I have performed supervised classification (Interactive supervised classification and Maximum Likelihood Classification) on landsat reflectance raster data to classify it into seven land use categories.

Now I want to know the area occupied by the mapped categories. Is there a function to do this in ArcGIS 10.1 Basic? The image is projected to UTM WGS 84 projection system. The area should be in sq km or sq m.


texts: Set a split policy/rule on a field without domain in ArcMap

texts:
Can you define a split policy or split rule on a field without a domain? When a domain is involved, you can set the default value of the field and chose a rule for when a feature is split or merged.

The setSplitRule exists, but the only two options are 'None' or 'Ratio'. 'None' sets the attributes to the copy of the original value.
http://resources.arcgis.com/en/help/main/10.2/index.html#//018z0000004v000000

I am looking to set a field to the Default Value (null) instead of copying the value from the original feature.


texts: How to get vector data from OSM to use in QGIS?

texts:
How do I get data from OpenStreetMap which I can then edit in QGIS layer by layer?

What is the best OSM output and what / how do I get to a suitable input for QGIS?
OR...is there another approach to this?


texts: Iterate - Raster Calculator - 8 Day MODIS LST to Monthly Avg

texts:
I have 8-day MODIS Land Surface Temperature data from 2009-2013. I have used Model Builder to extract the LST_day and LST_Night Subdatasets (in different model runs). What I am attempting to do next is take 3, 8-day TIFFS (combined, a month of data) and average them so it comes out to be 1 monthly avg. TIFF for each month from 2009 to 2013. I am unsure how I can select specific 8-day data to average in raster calculator using iterator - is there a way to specify or only select certain data? Is there a way to do this in ArcPy?

This is my current model:


texts: get external jsonp and convert coordinates in leafletjs

texts:
I need to solve two problems:

convert my coordinates from wgs84 to epsg28992
use external geojson from an jsonp url

Both problems are as far as I understand not supported by Leaflet itself. So I took a deep dive to found alternative methods.

The content of the jsonp looks like:

getmyjson({"type":"FeatureCollection","features":[
 { "type":"Feature",
   "properties":{"waarde":72.000000, "code":"002002", "klasse":1, "aant_min_stap":0, "stapgrootte":50.000000, "aant_stappen":4},
   "geometry":{"type":"Point", "coordinates":[129544.434902 , 512753.813614]}
 },
.....
 { "type":"Feature",
   "properties":{"waarde":54.000000, "code":"BDV032", "klasse":1},
   "geometry":{"type":"Point", "coordinates":[117763.299998 , 503647.849492]}
 }
]});

The first problem (convert coordinates) I solved by using jquery: Result

relevant code:

function getmyjson(data) {
      geojsonLayer.addData(data);
 }

 $.getJSON("rdtest.json", function(data) {
    var geojsonLayer = L.geoJson(data, {
      pointToLayer: function (feature, latlng) {
        latlngnw=RdWgs84(latlng.lng,latlng.lat)
    return L.circleMarker(latlngnw,{})
      }
    });
    geojsonLayer.addTo(map);
  });

The second problem (use external jsonp) I solved by using the plugin leaflet-ajax: Result

relevant code:

function getmyjson(data) {
      geojsonLayer.addData(data);
      geojsonLayer.addTo(map);
 }

 var geojsonLayer = new L.GeoJSON.AJAX(["wgs84test.jsonp"], {
      dataType:"jsonp",
      callbackParam:"getmyjson",
      });

What I want now is to combine both solutions and I have no clue how to do that.

Can anyone help me out ?

From ThomasG77 this is the solution to combine 1. and 2. Result

function getmyjson(data) {
   geojsonLayer.addData(data);
   geojsonLayer.addTo(mymap);

}

var geojsonLayer = new L.GeoJSON.AJAX(["rdtest.jsonp"], {
    dataType:"jsonp",
    callbackParam:"getmyjson",
    pointToLayer: function (feature, latlng) {
    latlngnw=RdWgs84(latlng.lng,latlng.lat)
    return L.circleMarker(latlngnw,{})
      }

});

But I will tweak further to look if I can use Proj4Leaflet to solve it.

Alternative for point 1 with use of Proj4Leaflet. Result

proj4.defs("EPSG:28992","+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs");
  $.getJSON("rdtest.json", function(data) {

    data.crs={
       "type": "name",
       "properties": {"name": "urn:ogc:def:crs:EPSG::28992"}
    };
    L.Proj.geoJson(data).addTo(map);
  });

Remaining question now is how I can combine the solutions of point 2 (jsonp) and the above.

Any suggestions ?

Final solution. Result

proj4.defs("EPSG:28992","+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs");

function getmyjson(data) {
   geojsonLayer.addData(data);
   geojsonLayer.addTo(mymap);
}

var geojsonLayer = new L.GeoJSON.AJAX(["rdtest.jsonp"], {
    dataType:"jsonp",
    callbackParam:"getmyjson",
    pointToLayer: function (data) {
          data.crs={
            "type": "name",
            "properties": {"name": "urn:ogc:def:crs:EPSG::28992"}
      };
    return L.Proj.geoJson(data)
        }
});


texts: Windows desktop application that allows to visualize, edit, download, upload and convert GPX files?

texts:
I am looking for a simple application that will allow me to visualize GPX data, edit and create new GPX tracks, areas and points and download/upload GPX data to and from devices. This application should have worldwide vector maps and imagery background coverage.

DNRGPS is good but cumbersome and lacks the ability to visualize and visually edit and create data.

Google Earth Pro does not export to GPX.