texts:
I have a WMTS service URL, which gives me three layers in projection - WGS84_Pseudo-Mercator (tileMatrixSet). It's very easy to add them using this guide -> https://openlayers.org/en/latest/examples/wmts.html . But I also need to add layer in "ETRS-TM35FIN" TileMatrixSet - as I understand  ETRS-TM35FIN is equal to EPSG:3067.

First solution that I found was to use raster reprojection from official docs http://openlayers.org/en/master/doc/tutorials/raster-reprojection.html but I DON'T want to change my View and set New View with another projection, because I have a lot of WFS layers with objects on the map and when I will change view to EPSG:3067 they all will disappear, because they live in EPSG:3857 system.

How could I load this layer in ETRS-TM35FIN (EPSG:3067) into my map with EPSG:3857 projection?

Please, anyone, help.


texts: While importing KML/KMZ rasters only colourful rectangles appear, although in right places

texts:
I'm trying to load KMZ rasters (file here) from Google Earth, but they are imported as mono-colour rectangles. They are in the right places though.

Global Mapper is suggested as a solution elsewhere (268166) but is not available for macOS.


texts: Does the store type Web Feature Server (NG) handle up to a specific number of features?

texts:
I am experiencing problems using the Store type: Web Feature Server (NG) in order to get Features from a remote GeoServer (2.5).

Although when initializing the store there is the possibility to choose:

  Maximum number of Features to retrieve (0 = no limit)

I made some tests by creating a layer with hundreds of thousand of features and it crashes. It seems like it works only with a few thousands of features but can not handle more.

GeoServer logs return this:

Could not acquire feature:java.net.SocketTimeoutException: Read timed out
java.util.NoSuchElementException: Could not acquire   feature:java.net.SocketTimeoutException: Read timed out

Caused by: java.net.SocketTimeoutException: Read timed out
org.geoserver.platform.ServiceException: Rendering process failed

Caused by: java.util.NoSuchElementException: Could not acquire feature:java.net.SocketTimeoutException: Read timed out

Caused by: java.net.SocketTimeoutException: Read timed out

In order to create the WFS Cascaded Server I followed the instructions here:

http://geoserver.geo-solutions.it/edu/en/adding_data/add_wfscascade.html

Leaving the following parameters with default values:

Connection and read timeout (ms): 3000
Feature buffer size: 10
Maximum number of Features to retrieve (0 = no limit): 0


texts: increment column value with an integer when it is a duplicate

texts:
I have a layer with around 1000 lines. Some Name are repeated several times and I wish to have every name unique. Thus I try in the field calculator to add an integer at the end of the name only when it is a duplicate, and of course I want to increment the integer value each time it is used.
As an example, in my layer I have in the Column 'Name':

Name_One
Name_One
Name_Two
Name_Two
Name_One

And I would like it to become

Name_One1
Name_One2
Name_Two1
Name_Two2
Name_One3

I tried first to identify the duplicate with the function

if(count("Name", "Name") = 1, 1, 0)

Then with the new column "uniqueness" that hold the value of the previous command I try

CASE
    WHEN 'uniqueness' = 0 THEN @a
    WHEN 'uniqueness' = 0 THEN @a = @a+1
    ELSE @a = 1
END

(@a is a project variable equal to 1)

But it doesn't work.

How can I do that?


texts:
I need to find how many routes are assigned to each account number in a feature class. I'm using search cursor twice, first to build a set of account numbers, second to get the number of times a route number is found associated with an account number. The second search cursor, I think, would be introduced after a for loop based on each value in the accounts set. As the search cursor is run through each value in the set, it should count the number of times a route number is found and then build a dict with the account number as the key and the number of different routes as the value.

Here is what I have so far:

import arcpy

# set workspace
arcpy.env.workspace = r"C:\RS_Data\Workspace\gisdb\layers.gdb"

# create input variable based on service location layer
inAccts = r"C:\RS_Data\Workspace\gisdb\layers.gdb\Accts"

# create empty set and dict
Accts = set()
AcctRts = {}

# Populate Accts with account numbers
with arcpy.da.SearchCursor(inAccts, "Account_ID") as cursor:
    for acct in cursor:
        Accts.add(acct[0])

# populate dict
for acct in Accts:
    with arcpy.da.SearchCursor(inAccts, "Exist_Route") as cursor:
        for row in cursor:
            RtCt = arcpy.GetCount_management(row[0])
            AcctRts[acct] = RtCt

for k in AcctRts.items():
    print k

And here is the error I'm getting:

Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 323, in RunScript
    debugger.run(codeObject, __main__.__dict__, start_stepping=0)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\__init__.py", line 60, in run
    _GetCurrentDebugger().run(cmd, globals,locals, start_stepping)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\debugger.py", line 654, in run
    exec cmd in globals, locals
  File "N:\Python\Completed scripts\MultiDay_Check_num_of_routes_per_Acct.py", line 1, in <module>
    import arcpy
RuntimeError: cannot open '051244-0001'
Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 323, in RunScript
    debugger.run(codeObject, __main__.__dict__, start_stepping=0)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\__init__.py", line 60, in run
    _GetCurrentDebugger().run(cmd, globals,locals, start_stepping)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\debugger.py", line 654, in run
    exec cmd in globals, locals
  File "N:\Python\Completed scripts\MultiDay_Check_num_of_routes_per_Acct.py", line 1, in <module>
    import arcpy
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\management.py", line 14483, in GetCount
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Rows: Dataset 021 does not exist or is not supported
Failed to execute (GetCount).


texts:
How to use the console.log() on openlayers3 objects?

For example if I log the map object in openlayers2 this is a result.
 Here i can observe the map, layers, baselayers etc. and its very helpful in development.

But with the new openlayers i.e. openlayers3 when i console log the map object i get something totally different in log. After googling a bit, i found this post, and tried logging the map object as map.getLayers() and i get this output
 which doesn't seems helpful compared to older version. So is there any way to log the openlayers objects so as to get the log clear as it used to in the version 2+. Any pointers related to this are most welcome.


texts: Traditional or Vintage mapping theme in ArcMap

texts: Is there a labeling or theme extension that can be downloaded for ArcMap that would help to make my maps look Vintage like historic maps? I've looked at some of the ESRI extensions and basemaps online, but there are so many I wouldn't really know where to start. I was hoping to zero in on one that someone had utilized and would recommend. Any suggestions?


texts: Is there a way for me to load a KML files locally in Openlayers

texts:
Recently I posted this question KML file seems to load in Openlayers without error but doesn't show vectors because I was having trouble getting KML vector data to show up. After searching around a bit it seems there is an issue with loading KML files locally.

I am developing a stand alone Application that will run on the desktop or on IOS. I am running openlayers via a WebKit view. I am loading my openlayers page from the local file system witch is part of the applications files when it's installed.

My guess is there is some sort of issue loading local files. Is there a way around this? I would even consider converting all my KML files if needed. Ideally I would like to be able to point to either a local file or a remote kml file and load it.

Here is the code I am using. It's from the Openlayers kml example that's online. The only difference is I am trying to load locally without running the html page on any sort of server.

var map = new OpenLayers.Map({
    div: "map",
    layers: [
        new OpenLayers.Layer.WMS(
            "WMS", "http://vmap0.tiles.osgeo.org/wms/vmap0",
            {layers: "basic"}
        ),
        new OpenLayers.Layer.Vector("KML", {
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.HTTP({
                url: "kml-files/lines.kml",
                format: new OpenLayers.Format.KML({
                    extractStyles: true,
                    extractAttributes: true,
                    maxDepth: 2
                })
            })
        })
    ],
    center: new OpenLayers.LonLat(-112.169, 36.099),
    zoom: 11
});


texts: ArcGIS Server directory use to store intermediate files

texts:
I have a tool which does the following work:

Connects and fetches data from hdf using impala
Create intermediate csv file.
from the csv file it creates custom feature classes.

Currently when running it as a standalone tool I have used absolute paths in the local disks and it works fine. I have used random number to generate different work directories to store the data.

Now I need to publish this tool as a geoprocessing service. So I need to access the server directories to save my intermediate files.

As I read, each time a tool is run it creates a folder with a unique id within arcgisjob folder. If I understood right, I need access to that directory.

Please suggest how to use the folder.

I am using ArcGIS 10.1 Server for Windows.


texts: Zoomable world map with small countries highlighted

texts:
I'd like to create a zoomable world map where countries/territories are color-coded.  This is a pretty standard request and Carto has a nice example.
However, in standard maps like this, small countries are practically invisible unless you zoom in, so I'd like to highlight them the way Wikipedia often does:

Larger version here, and note that small countries have a circle superimposed on top.
Is there an existing library/package/... that can do this in a popular tool like Carto?  Bonus points if the circles go away when you zoom in close enough.  The closest I could find with some searching was Using 'Style by Value' with Point Geometries, but this uses dots instead of country polygons, which is not what I'm looking to do.


texts: Using Vector Tile Maps Offline in a Hybrid App (ArcGIS API for JavaScript)

texts:
I will create a hybrid map by using ArcGIS API for Javascript. I want to use my own hosted vector tile basemap as offline. What do you advise?

Are there any ways to use ESRI Vector Tile Package (.vtpk) file as offline without serving on portal or Arcgisonline?


texts: How to call Exelis ENVI tools with python and ArcGIS?

texts:
I am writing a script for geopocessing landsat data using ArcGIS 10.3. But ArcGIS doesn't have tools for landsat calibration (converting from DN to reflectance) and ENVI have this tool called "Radiometric calibration".

There is a way to call Exelis ENVI 5.2 tools using python?


texts: Using an ArcGIS Model Builder iterator to merge all shapefiles into one layer?

texts:
I have a GeoDatabase with a handful of shapefiles. The contents are likely to change often, so I want to build a tool that allows me to merge all the shapefiles in the GDB into one layer.

This is what I currently have:

However, I am not expecting this to work, as it seems like it would just grab one file each iteration and 'merge' it with itself. Is there a workaround so that it grabs all files from the GDB to merge them? I am using ArcGIS Desktop 10.5.


texts: Editing CharacterSpacing annotation attribute in ArcMap?

texts:
I want to change the character spacing for a single (selected) feature/record in a file geodatabase annotation feature class.

I understand that the CharacterSpacing attribute controls that, but changing the default (0) to other values (positive or negative) does nothing. By comparison, editing CharacterWidth does work as expected.

Any ideas?

I am running ArcGIS Desktop 10.5 Advanced on Windows 10.


texts: Convert lines to polygons with shapely polygonize

texts:
I have a line shapefile that I am trying to convert to polygons with python using fiona and shapely.

inFile = '/content/drive/My Drive/polyline.shp'
outFile = '/content/drive/My Drive/polygon.shp'
polylines = fiona.open(inFile)

I know i can loop though the data to look at it

for row in polylines:
    print  (row)

and looking at the data

{'type': 'Feature', 'id': '1', 'properties': OrderedDict([('ID', 1), ('Feld1', 189.02), ('Feld2', 'Kürbis')]), 'geometry': {'type': 'LineString', 'coordinates': [...........]}}

{'type': 'Feature', 'id': '2', 'properties': OrderedDict([('ID', 2), ('Feld1', 200.87), ('Feld2', 'Aprikose')]), 'geometry': {'type': 'LineString', 'coordinates': [...........]}}

I am wondering how I can use shapely polygonize with this data?


texts: Grid parameter restriction of mpl_toolkits.basemap.maskoceans

texts:
I have plotted the image below with the use of matplotlib and mpl_toolkits.basemap. However I am not satisfied with the result. The image is still very blocky near the coast line. The grid parameter of the maskoceans method seems to be the issue, however the API doesn't allow a lower grid size. How can I work around this restriction and why did they made this restriction? I guess you don't have this problem for larger countries.

Here a part of the code I used:

xArray = np.array(xArray)
yArray = np.array(yArray)
heightArray = np.array(heightArray)
xmin, xmax = 3.3, 7.3
ymin, ymax = 50.7, 53.5

numIndexes = 1000
xi = np.linspace(3.3, 7.3, numIndexes)
yi = np.linspace(50.7, 53.5, numIndexes)
zi = ml.griddata(xArray, yArray, heightArray, xi, yi, interp='linear')

m = Basemap(projection = 'merc', llcrnrlon = xmin, llcrnrlat = ymin, urcrnrlon = xmax, urcrnrlat = ymax, resolution='h')
m.drawcountries(linewidth=1)
m.drawcoastlines(linewidth=1)

xx, yy = m(np.array(xi), np.array(yi))
xi, yi = np.meshgrid(xi, yi)
zi = maskoceans(xi, yi, zi, resolution = 'h', grid = 1.25, inlands=True)
m.drawlsmask(land_color='grey',ocean_color='white', grid=1.25, resolution='h')
m.pcolormesh(xx, yy, zi, cmap = plt.get_cmap('jet'), vmin=0, vmax=100)
plt.colorbar()

To clarify what I would like even more I added the image below. This image was edited by a photo editor and is not something what matplotlib produced directly.


texts: I am looking for a way to add Google Earth imagery to ArcMap 10.  Over a year ago, a GIS.SE member mentioned Arc2Earth was the only way to accomplish this (04/2011 thread).  Another somewhat convoluted approach is described here.  Has there been any progress in bringing Google Earth imagery to ArcGIS?


texts: Predict a pair of geographical coordinates from their equivalent in flat coordinates in Google Colab

texts:
I try to predict unknown geographic coordinates from meters coordinate in an excel table. (x and y). something like this:

import pandas as pd
import numpy as np
from sklearn import linear_model
from sklearn.metrics import mean_squared_error, r2_score
from math import radians
import matplotlib.pyplot as plt

First I will connect to drive server

from pydrive.auth import GoogleAuth
from pydrive.drive import GoogleDrive
from google.colab import auth
from oauth2client.client import GoogleCredentials
# Authenticate and create the PyDrive client.
auth.authenticate_user()
gauth = GoogleAuth()
gauth.credentials = GoogleCredentials.get_application_default()
drive = GoogleDrive(gauth)

Them, I get the data

import pandas as pd
link1 = 'https://drive.google.com/open?id=1rOESJV89SzIDiQku88MsIXA6_WcyB0Ol'
wdir, id = link1.split('=')
id

downloaded = drive.CreateFile({'id':id})
downloaded.GetContentFile('data.xlsx')
print('The files was downloaded')

df = pd.read_excel('data.xlsx', index_col=0)
df.columns = ['x','y','W','N']
df['Id'] = ['n_{}'.format(i) for i in range(len(df))]
print(df.info())
df.head()

As you can see, I get missing data in W and N columns, that there is the geographic coordinate:

However, the columns x and y are complete. Is it possible to predict the N and W columns use x and y columns? I was using the next pace of code (lineal regression) for my propose:

sample = df[df.N.notnull() & df.W.notnull()].copy()

regr = linear_model.LinearRegression()

y = sample[['W','N']].values

X = sample[['x','y']].values

regr.fit(X,y)

r2_score(sample[['W','N']].values, regr.predict(sample[['x','y']].values))

If you gonna see my notebook this is the link:
Notebook in Google Colaboratory


texts: Is it possible to play records which has trajectories and time stamps as live traffics. I did investigate the World Wind SDK but I can not see special classes for aviation or Tracks and trajectory calculations.


texts: Is there a public CanTopo Web Map Tile Service (WMS / WMTS) available?

texts:
Does anyone know of a public CanTopo WMS / WMTS available anywhere? I'm assuming CanTopo is more modern / up-to-date than Toporama?

Alternatively, what would be your suggestions for the most pain free and cheapest version to host such a service myself?

The data seems to be freely available here:
http://www.nrcan.gc.ca/earth-sciences/geography/topographic-information/data/11042


texts: Issues with ArcGIS Pro and Python for Labeling Point Shapefiles

texts:
I am trying to take two separate columns of data and create a stacked label where the bottom line (Benzen_1Q_2019) will be red and a larger font, while the top line (Well_ID) will be a dark grey color and a smaller font. I have symbolized the labels already and am now working on using Python to change the color and the font size. I have gotten the majority of it done, but the one thing I am struggling with is the column (text field) I am using for the second line (Benzen_1Q_2019) has values that are either a whole number or (-) or <1. When I type out out the Python expression, it works correctly for the values that are either (-) or a whole number, but it does not work for the values that are <1.

Any tips for how to make this work correctly?

My expression:

[Well_ID] + '\n' + "<CLR red ='255'><FNT size ='7'>" + [Benzen_1Q_2019] + "</FNT></CLR>"

Attribute Table
Python code and labels]2


texts:
I need all labels to be placed outside of the boundary (polyline features).  Most of the labels are placed outside the boundary line, while others (marked with red color) are placed inside the boundary.

What is the best way place all labels on the outside of the polyline?


texts: How do I use python to batch process the reprojection of a group of shapefiles?

texts: My client would like his shapefiles reprojected and I would like to get a script that can batch process all the files located in the same folder. I am new to Python and ArcPy but am very interested. We are using ArcMap10 as our platform software. Any help or advice would be greatly appreciated.


texts: Error 000725 importing Python (Arcpy) script into tool of standard toolbox?

texts:
I have created a Python script. When I run the script as a standalone by make this a hard set path instead of sys.argv1, the script really runs fine.

When I import this script as a python tool box then I get the error at output workspace that dataset is already existed. But there is no dataset in the GDB. I have attached the error file and the script file as well. I think i have set up the parameters  right as well.


texts:
I'm working on code, and I'm trying to select features by location. After the features are selected, I want the cursor will add "Yes" to "NearPoint3" field, only to the features that were selected.
Right now the word "Yes" is added to all features.
The relevant part is on the last section.
import arcpy
import os

input_folder = r'C:\Users\roniy\Documents\UN\arcpy_final\final_script\Early_childhood_frameworks.shp'
OutPut_Folder = r'C:\Users\roniy\Documents\UN\arcpy_final\final_script\out'
OutPut_intersect = r'C:\Users\roniy\Documents\UN\arcpy_final\final_script\out_inter'
out_copy = r'C:\Users\roniy\Documents\UN\arcpy_final\final_script\out_copy'
OutLayerSelection = r'C:\Users\roniy\Documents\UN\arcpy_final\final_script\sell_fin'
arcpy.env.workspace = input_folder
fc = "Early_childhood_frameworks"
os.path.join(OutPut_Folder, fc)
Buffer = arcpy.Buffer_analysis(fc, OutPut_Folder, "100 Feet")
Intersect = arcpy.Intersect_analysis(Buffer, OutPut_intersect)

arcpy.MakeFeatureLayer_management(Buffer, "buffer")
arcpy.MakeFeatureLayer_management(Intersect, "intersect")
selection = arcpy.SelectLayerByLocation_management("buffer", "INTERSECT", "intersect")
LayerSelection = arcpy.CopyFeatures_management("buffer", out_copy)
clip = arcpy.Clip_analysis(fc, LayerSelection, OutLayerSelection)

#------------------------------------------------------

NewField = arcpy.AddField_management(fc, "NearPoint3", "TEXT")
arcpy.MakeFeatureLayer_management(clip, "cliped_p")
arcpy.MakeFeatureLayer_management(fc, "layer")
selection2 = arcpy.SelectLayerByLocation_management("layer", "INTERSECT", "cliped_p")
cursor = arcpy.da.UpdateCursor(fc, ["NearPoint3"])
for row in cursor:
    cursor.updateRow(["Yes"])


texts: Getting multiple dynamic custom markers on map using MapBox in MySQL and PHP?

texts:
Can anyone help me to get multiple custom markers on mapbox using lat/long from the database MySQL?

I want to use 2 types of markers on the map, for which I have two separate table of lat/longs in the database.

I already done it by Google map but now I would like to do it with MapBox. My current code to show those markers in a map:

  Google Maps Multiple Markers

  <?php
    $a= array();
    $i=0;
    $sql= mysqli_query($con,"SELECT * FROM mapdata");

    while( $row =mysqli_fetch_array($sql)){
      $b=array();

      $b[]=$row["address"];
      $b[]=$row["latitude"];
      $b[]=$row["longitude"];
      $b[]=$row["id"];
      $b[]=$row["name"];
         $a[$i]=$b;
         $i++;
     }
  ?>

  <?php

    $z = array();
    $n=0;
    $sql2= mysqli_query($con,"SELECT * FROM mapdatablue");

    while( $row2 =mysqli_fetch_array($sql2)){
      $y=array();

      $y[]=$row2["address"];
      $y[]=$row2["latitude"];
      $y[]=$row2["longitude"];
      $y[]=$row2["id"];
      $y[]=$row2["name"];
      $z[$n]=$y;
      $n++;
     }

   ?>

<script type="text/javascript">
   var locations = <?php echo json_encode($a);?>;
   var locationsb = <?php echo json_encode($z);?>;
   var map = new google.maps.Map(document.getElementById('map'), {
       zoom: 6,
       center: new google.maps.LatLng(22.00, 96.00),
       mapTypeId: google.maps.MapTypeId.ROADMAP
   });

   var infowindow = new google.maps.InfoWindow();
   var image = '/marker/map.png';
   var image2 = '/marker/Map1.png';

   var marker, i;

   for (i = 0; i < locations.length; i++) {
      marker = new google.maps.Marker({
      position: new google.maps.LatLng(locations[i][1], locations[i][2]),
      offset: '0',
      icon: image,
      title: locations[i][4],
      map: map
   });

   google.maps.event.addListener(marker, 'click', (function(marker, i) {
      return function() {
        infowindow.setContent(locations[i][0]);
        infowindow.open(map, marker);
        }
      })(marker, i));
   }

   var marker1, n;
      for (n = 0; n < locationsb.length; n++) {
         marker1 = new google.maps.Marker({
           position: new google.maps.LatLng(locationsb[n][1], locationsb[n][2]),
           offset: '0',
           icon: image2,
           title: locationsb[n][4],
           map: map
        });

     google.maps.event.addListener(marker1, 'click', (function(marker1, n)
     {
       return function() {
          infowindow.setContent(locationsb[n][0]);
          infowindow.open(map, marker1);
       }
     })(marker1, n));
   }

  </script>
 </body>
</html>

**Output of this code **


texts: Merge Featureclasses if there is more than one in the FeatureClassSet

texts:
i want to check through all my gdb if a FeatureClassSet contains more than one FeatureClass. If so merge them and save the result in an other gdb otherwise just copy them to the gdb. It sounds really simple to me - but it is not working. So far I've got this:

import arcpy, os
from arcpy import env

workspace = arcpy.GetParameter(0)
EK = arcpy.GetParameter(1)

outWorkspace = r"C:\...\20_Ges.gdb"
Zieldataset = os.path.join(outWorkspace, EK)

fc = arcpy.ListFeatureClasses(feature_dataset=EK)
for item in fc:
    item.count(<1):
        arcpy.FeatureClassToFeatureClass_conversion (workspace, Zieldataset)
    else:
        arcpy.Merge_management(fc, Zieldataset)

I connected the script to a model that iterates the variable EK and sets the workspace. The error in the script is in the line that starts with "item.count".


texts: Create list of [X, Y] coordinates within field calculator

texts:
Lately i've been working on a script that must calculate a value in ArcGIS field calculator based on all vertex coordinates of a polyline. But i am struggling to create a list of a polylines coordinates whitin the field calculator. The list is not meant to fill the field values, but it make other calculations possible.

This is how the list must be like:

coordinates = [[534359.1787999999, 6944339.3809], [534298.1876999997, 6944310.3992], [534282.0721000005, 6944304.4425]]

I have found some material about it on ArcGIS Resource Center, it works when i run their script through the command line and append to a list, but it doesn't work at all in field calculator. I've tried to write a code for the field calculator but it didn't work either, it look like this:

def listCoordinates(shape):
    coordinatesList = []
    n = 0
    for pnt in shape.getPart(n):
        xn = pnt.X
        yn = pnt.Y
        coordinatesList.append([xn , yn])
        n += 1
    return coordinatesList

and the expression (this function will be inside another function that return a unique value):

listCoordinates(!shape!)

Is this possible? Am i missing something? Any tips will be welcome! :)

It turns out that arcpy.getPart() do not return an array in field calculator. So you can use arcpy.Array() to get to the array. I was able to get what i wanted with:

def getCoords(shp):
    coords = []
    n = 0
    for part in arcpy.Array(shp.getPart(n)):
        for p in part:
            coords.append([p.X , p.Y])
        n += 1
    return coords


texts: Openlayers Measure Line and Area at the same time from the same button

texts:
I have been using separate buttons to measure lines and areas in my Openlayers maps (as per the Openlayers examples) but would ideally like to get the output for both length and area at the same time. The ideal solution would be to have a solid line for the line measurements and a dashed line connecting the start and end point to symbolize the polygon. Then simply output the two answers to a div on the website. As a bonus, store each line length in an array that can be output to a menu as well.

My question is two fold, is it possible, and does anyone have / know of an example? I've searched long and hard and can't find an attractive solution to this.

Thanks for your help,


texts:
I just started using gdal. I'm trying to get longitude and latitudes from pictures. This is what I have so far. This gives me an attribute error. Anyone know how to fix the attribute error?
    from osgeo import osr
    import numpy
    import gdal
    ds = gdal.Open(r"path")
    prj = ds.GetProjectionRef()
    oldcs = osr.SpatialReference()
    oldcs.ImportfromWkt(prj)

Error:
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-1-b32399e41073> in <module>
      7 prj = ds.GetProjectionRef()
      8 oldcs = osr.SpatialReference()
----> 9 oldcs.ImportfromWkt(prj)

D:\Anaconda\lib\site-packages\osgeo\osr.py in <lambda>(self, name)
    354     __setattr__ = lambda self, name, value: _swig_setattr(self, SpatialReference, name, value)
    355     __swig_getmethods__ = {}
--> 356     __getattr__ = lambda self, name: _swig_getattr(self, SpatialReference, name)
    357     __repr__ = _swig_repr
    358

D:\Anaconda\lib\site-packages\osgeo\osr.py in _swig_getattr(self, class_type, name)
     78     if method:
     79         return method(self)
---> 80     raise AttributeError("'%s' object has no attribute '%s'" % (class_type.__name__, name))
     81
     82

AttributeError: 'SpatialReference' object has no attribute 'ImportfromWkt'


texts:
Using the new (in 2.2) 'Fill Ring' tool and it doesn't seem to work correctly in all cases. It seems that if I try and digitise the new ring close to the polygon boundary, TWO (or more) geometries will be created (and will get prompted to fill the attributes for), one for each adjacent polygon.

If I digitize the ring around the polygon center everything works fine.

Same behavior for shapefiles and Postgres layers that I tried.

The screenshot below shows an example. The attributes form displayed is the SECOND one and refers to the attributes of the polygon on the LEFT of the polygon I added the ring to.

Anyone got any ideas? Is this a bug?


texts: I did a bunch of "housekeeping" recently—deleting old unneeded attribute fields, renaming other fields, etc.. Now I find that all of the labels that I had previously manually repositioned  seem to have shifted down a small amount (by an amount approximately equal to the size of the label letters). Any idea what might have caused this, and how I can correct it (other than by manually readjusting all of them (hundreds))? And more importantly, so I can avoid causing this to happen again? (Using QGIS 3.4)


texts: Why does GeoServer/OpenLayers show only points from my PostGIS?

texts:
I am reading the PostGIS in Action book now and I have a problem with one example:

I have created a table:

CREATE TABLE my_geometries (id serial NOT NULL PRIMARY KEY, name varchar(20));
I have added some points to this table:

SELECT AddGeometryColumn('public','my_geometries',
'my_points',-1,'POINT',2);

INSERT INTO my_geometries (name,my_points)
VALUES ('Home',ST_GeomFromText('POINT(0 0)'));

INSERT INTO my_geometries (name,my_points)
VALUES ('Pizza 1',ST_GeomFromText('POINT(1 1)')) ;
And some linestrings as well:

SELECT AddGeometryColumn ('public','my_geometries',
'my_linestrings',-1,'LINESTRING',2);

INSERT INTO my_geometries (name,my_linestrings)
VALUES ('Linestring Open',
ST_GeomFromText('LINESTRING(0 0,1 1,1 -1)'));

INSERT INTO my_geometries (name,my_linestrings)
VALUES ('Linestring Closed',
ST_GeomFromText('LINESTRING(0 0,1 1,1 -1, 0 0)'));

The problem is that for some reason I do not see the linestrings in the OpenLayers window of GeoServer, but all the points are on their place.

The GeoServer's 'Edit layer' page for my_workspace:my_geometries shows Feature Type Details having both my_points and my_linestrings columns.

What is the problem here? How can I see both points and linestrings that I have added to my database in GeoServer's OpenLayers window?


texts: Create multivariable NetCDF from several Raster files

texts:
Is it possible to create a NetCDF with several variables from a collection of rasters (same spatial bounds)?

I've tried the Arcpy tool RasterToNetCDF, but it seems thats it`s only able to transform one raster each time.


texts:
Is it possible to divide a QGIS scale bar into different sized segments as you can in MapInfo (as with example given).

Or is it possible to select different style of scale bars in QGIS?

If it's not already obvious I'm normally a MapInfo user.


texts: How to search for visible layers in a map?

texts:
I'm using ArcGIS runtime SDK for .Net, I have an application with only a

  ArcGISLocalDynamicMapServiceLayer

which contains a number of feature classes.

How can I get all the layers included in the Dynamic map service layer IN CODE in order to change the visibility of some of the layers individually?

XAML code:

<esri:Map x:Name="MyMap" UseAcceleratedDisplay="True" Cursor="Cross" Extent="330671.7091,2854689.5845,344178.0091,2862125.3845" MouseClick="QueryPoint_MouseClick" IsEnabled="False" IsLogoVisible="False">

            <esri:ArcGISLocalDynamicMapServiceLayer ID="MYLAYERS"   Path="C:\Users\Ahmad Elsayed\Desktop\test app\MODREN UI APP\FinalMap.mpk" />

  </esri:Map>


texts: Find points within polygon with PostGIS really slow

texts:
From some millions of points I am trying to select those points that are within a polygon and save them into a new table. Both datasets with spatial index and both on projected coordinates.

The query I am trying with is:

CREATE TABLE result AS
SELECT t1.point_id, t1.geom
FROM points t1
INNER JOIN boundaries t2
ON ST_WITHIN(t1.geom, t2.geom)
WHERE t2.country LIKE "ENGLAND"

The query has been running for hours.

On the same machine I have done exactly the same thing using ArcMap in 25 minutes.
I really don't understand this. I was under the impression that doing geoprocessing tasks directly on a RDBS (PostgreSQL or SQL Server) was the fastest option (provided the query is the right one).

Is there anything wrong with the query? Is there any way to optimise it? Am I missing something here?

PS: I tried with ST_Intersect instead of ST_WITHIN with the same result.

As suggested in 'Comments', I have tried with FROM being the smaller table (table with 3 polygons).

CREATE TABLE result AS
SELECT t1.point_id, t1.geom
FROM boundaries t2
INNER JOIN points t1
ON ST_WITHIN(t1.geom, t2.geom)
WHERE t2.country LIKE "ENGLAND"

and also

CREATE TABLE result AS
SELECT t1.point_id, t1.geom
FROM boundaries t2
INNER JOIN points t1
ON ST_Contains(t2.geom, t1.geom)
WHERE t2.country LIKE "ENGLAND"

Unfortunately, these two options seem to be pretty slow too. 2 hours running and still processing. I must admit that the polygon geometry is pretty complex (England boundary). But still, as indicated in the descriptions, other GIS tools take 25 minutes in doing the same job.

Explain analysis: (from last query written above, the one with "ST_Contains")

"Nested Loop  (cost=111.22..26768.81 rows=3007 width=44)"
"  ->  Seq Scan on gb_bdry t2  (cost=0.00..17.88 rows=3 width=32)"
"        Filter: ((country)::text ~~ 'ENGLAND'::text)"
"  ->  Bitmap Heap Scan on final_def_sgm_centroids t1  (cost=111.22..8906.96 rows=1002 width=44)"
"        Recheck Cond: (t2.geom ~ geom)"
"        Filter: _st_contains(t2.geom, geom)"
"        ->  Bitmap Index Scan on final_def_sgm_centroids_geom_idx  (cost=0.00..110.96 rows=3007 width=0)"
"              Index Cond: (t2.geom ~ geom)"

Different approach following suggestion from losbaltica:

First we limit the number of points with &&:

CREATE TABLE temp1 AS
SELECT t1.point_id, t1.geom
FROM points t1
WHERE t1.geom && (SELECT geom FROM boundaries WHERE country = 'ENGLAND')

Then we do the actual intersection:

CREATE TABLE temp2 AS
SELECT t1.point_id, t1.geom
FROM points t1,
    (SELECT geom FROM boundaries WHERE country = 'ENGLAND') AS t2
WHERE ST_Intersects(t1.geom, t2.geom)


texts: Stream network doesn't match the reality?

texts:
I am trying to create flow accumulation grid in ArcMap using Arc Hydro.
I have the DEM (25mX25m) and the stream network shapefile.

I have done the following

1)Resampled the DEM to 5m;

2)Burnt the stream network in the DEM using DEM reconditioning;

3)Filled the sinks;

4)Created flow direction grid; and

5)Created flow accumulation grid.

The flow accumulation grid I have created was matching well with the reality in most of the catchment. However, in the middle part of the catchment, two streams are connected while they are not connected in reality. In addition, the catchment outlet created by Arc Hydro was not matching with the reality as illustrated below.

Any recommendations how can I make the flow accumulation grid created by Arc Hydro match the reality would be appreciated?

Thanks in advance for your help and time.


texts:
I've got two GIS tables, let say A and B. I would like to do the difference of A using B. While it's pretty straightforward to do it with QGIS (using geoprocessing->difference), I've got much more troubles to do it with PostGis.

Below a visual example :

1/ My two layers (in yellow and orange)

2/ The difference (in pink) between my yellow layer and my orange layer (using QGIS geoprocessing->difference). Only the part of orange that does not intersect my yellow layer are kept

This is my SQL query :

SELECT  A.id,
    ST_Union(
        COALESCE(
            ST_Difference(
                A.the_geom,
                B.the_geom
            ),
            A.the_geom
        )
    ) As the_geom
FROM A
LEFT JOIN b ON ST_Intersects(A.the_geom, B.the_geom)
GROUP BY A.id;

This approach almost work, the only limit is when an object A intersects more than one B object. In this case, the difference function does not work at all.

As you can see in the picture below (difference in purple) the big shape in the middle remains the same instead of being cut by the two yellow geometries.

So my idea was to do a geometry union of the geometry from B :

SELECT  A.id,
    ST_Union(
        COALESCE(
            ST_Difference(
                A.the_geom,
                ST_Union(B.the_geom)
            ),
            A.the_geom
        )
    ) As the_geom
FROM A
LEFT JOIN B ON ST_Intersects(A.the_geom, B.the_geom)
GROUP BY A.id;

But Postgresql complains about the fact that aggregate functions cannot be nested.

Do you have any idea how I could do the same QGIS difference process but with PostGIS?


texts: Moving a Map Service Group Layer out of a Group Layer with ArcPy

texts:
I am struggling to move web services (tiled basemaps) out of group layers into the data frame. I've used arcpy.mapping.ListLayers to compose a list of layers, then I've ascertained what is a group layer with .isGroupLayer and I've tested the group layers name and then what it "supports". However I failed miserably to pick up weird web service group layer and move them. I've tested arcpy.mapping.AddLayer among a handful of things.

Any ideas are welcome. My desired outcome is to have no web services in group layers.

import arcpy
mxd ="C:\\Data\\native.mxd"
OriginMxd = arcpy.mapping.MapDocument (mxd)
df = arcpy.mapping.ListDataFrames(OriginMxd)[0]

for OPFSLyr in arcpy.mapping.ListLayers(OriginMxd, "", df):
    if str(OPFSLyr) == "Basemap": # if  it's named basemap
        if OPFSLyr.isGroupLayer:
            for subLayer in OPFSLyr:
                print subLayer
                arcpy.mapping.AddLayer(df, subLayer, "BOTTOM")
                OriginMxd.save()


texts: TPK creation failed - ERROR 001117: Map description property must be set Failed to execute

texts:
I'm using ArcMap10.6 to make a tile package. I've filled out the description and map title in the map document properties but I get the this error:

ERROR 001117: Map description property must be set Failed to execute

Which says the description in the map document properties must be populated:

Any clue why this would be happening?


texts: Seeking vector data for wind speed and operating turbines in UK?

texts:
I have been familiarising myself with QGIS and the free information Ordnance Survey provide. Though I'm struggling to find further data specific to my needs.

I want to produce a map plotting renewable energy in the UK (or at least operating wind turbines) against locations of operating data centres, preferably on a background of wind speeds.

Are there places where I access this data?

By so far I have seen only these sources:

This website provides locations of data centres: https://www.datacentermap.com/
This renewable energy: https://ukdataexplorer.com/renewables/
This has wind speed but when I import it is all the same colour: https://www.renewables-atlas.info/


texts:
I have a raster that I clipped to a polygon using the image analysis window (regular raster clip doesn't do what I want) but when I go to save the temporary raster to a hard file I get the same results as using the regular raster clip tool.  Is there any way to save the temporary image analysis clipped raster so it appears exactly the same way in a hard raster file?

Using the raster clip tool:

Temporary file produced when using image analysis (desired result when saving as a hard raster file):

When saving temp image analysis to a hard file:


texts: How to update attribute of selected features in ArcMap with a shortcut?

texts:
I have to edit/move/create/delete a no of polygons based on their respective comparison with reality. For all the polygons that I've moved, I have to update a field called remark as "U". This is really annoying to switch everytime to the attribute window and update the same.

Is there any shortcut/add-in/extension that could enable this step in a single click(shortcut key) or automated if possible?

Liscence is Basic for ArcMap 10.3


texts: Is there a style set out there that has state shields for labeling highways? For example: one of the labels has the Interstate Highway symbol, but is there a style set that is state specific, such as FM roads or Texas State Highway roads?


texts: CARTO.js / ReactJS: Showing and hiding layers on toggle

texts:
I have a the following layer:

class DataLayer extends Component {
  static contextTypes = {
    map: PropTypes.object, // Leaflet map
  };

  componentDidMount() {
    const { client } = this.props;

    // Craete source, styles and layer with the given props
    const cartoSource = new carto.source.Dataset('dataset');
    const cartoCSS = new carto.style.CartoCSS(style);
    const layer = new carto.layer.Layer(cartoSource, cartoCSS);

    // Add them to the client and to the map
    client.addLayer(layer);
    client.getLeafletLayer().addTo(this.context.map);
  }

  // ... missing methods to handle styles/source updates

  render() {
    return null; // No need to render anything :)
  }
}

export default DataLayer;

Then, I have the following component that shows the Layer component:

class Canvas extends Component {
  state = {
    center: [40.6552398976147, -73.959729361534119],
    zoom: 11,
    checked: false
  };

  constructor(props) {
    super(props);

    // Setup the client in the contructor with our user and apiKey
    this.cartoClient = new carto.Client({
      apiKey: "APIKEY",
      username: "LOGIN"
    });
  }

  onClick(e) {
    this.setState({
      checked: !this.state.checked
    });

  }

  render() {
    const { center, zoom } = this.state;

    return (

                  <input
                    type="checkbox"
                    onClick={this.onCheck.bind(this)}
                    checked={this.state.checked}
                  />

          <Map center={center} zoom={zoom}>

            <Basemap attribution="" url={CARTO_BASEMAP}/>

            {this.state.checked ? <DataLayer client={this.cartoClient}/> : null}
          </Map>

    )
}
export default Canvas;

Everything works great for the first try, when the checkbox triggers checked to be true and the DataLayer component mounts. However, when I uncheck the checkbox and expect the DataLayer component to be unmount, it doesn't. I believe the reason is a missing method to handle it in the DataLayer component once it's been updated.

How to make DataLayer component to "unmount" when the toggle is set to false?


texts: Converting multiple shapefiles into raster using ArcGIS ModelBuilder?

texts:
I want to convert multiple shape files into raster files using iteration in ModelBuilder in ArcGIS. Each time I am facing the same error, which is showing:

  Error executing function. A column was specified that does not exist.

My model is attached.


texts: How can I create a map style for OSM data in QGIS 3

texts:
I am new to QGIS (Version 3.2.2) and basically what I want to archive is to load OSM data and apply my own map style to it. However once I have defined my style I would also like to be able to load new OSM data and instantly select the same features (roads, rails, etc.) and apply my custom style.

The workflow I have come up so far is:

load the data according to https://gis.stackexchange.com/a/101746/120164
creating my style by creating layers, which I do as follows:

connect to the database
expand the database in the browser
double click on the type of features I want in my new layer (points, lines, polygons, etc.)
add a filter to the new layer, so only features with the key/tag I want are in the layer

From a software developers point of view, I would now just switch the database connection and the queries could be applied to another database and so a map of another region would be rendered without having to do much.
However it seems that the layers are somehow linked to the database.

Is there an easy way to just take the style/selection of features and apply it to other mapdata?


texts: Convert Image Band Values from Array into Number

texts:
I am trying to convert the pixel values that appeared to be an array (see the snapshot below) into numeric. In its current state, every time I tried implementing reduction, I got an error:

Input to Reducer.mean must be a numeric scalar, not
Type<Float<dimensions=1>>.

Here's a snapshot of the console when I printed the data:


texts: Center generated atlas map to features of different layer

texts:
I got a problem with setting an Atlas. I have a polygon layer with around 100 polygons, and a line layer.

What i need to do is to generate a map for each polygon and show on the fixed scale (1:500) those specified lines.

I configured atlas to generate those maps based on polygon layer, but when polygon is large, lines are not visible because atlas centering map to the centre of polygon. Is there a way to generate map centered to those lines?

Im using a "inverted polygons" as a style of Polygon layer becouse im showing only one specific polygon that im reffering to on atlas using "$id = @atlas_featureid" in expression.

Got some screenshots to show what is in my mind:

Screen 1 - Nicely generated map

Screen 2 - Map generated in 1:1000 scale showing the problem and the boundary i would like to show while in 1:500

Screen 3 - Badly generated map


texts: Nesting styles for increased drawing speed

texts:
I have zoning plans for a municipalities ( for example 5700 features in this one) and I have created a QGIS-style to make it look like the offical style.

There are 72 colours/patterns for the different zones (coarse version), and  at least twice the number of zones in a more detailed style. Some of the zoning codes share the same patterns/colours as other codes.

All the codes are four digit so I tried to regex the zoning-code and categorised them like shown on the image below.

Running this style on this number of features takes some time every time I do anything in QGIS.

So the question i: I was wondering if nesting the zone-styles would help with performance, or if this would just end up using more (simple check vs regex for the category -> sub-category)


texts: Feature id generation in spatial database

texts:
I am using a view

SELECT a.id,
    c.nombre AS "UEB",
    b.nombre AS "Unidad",
    a.geom
   FROM genetica.multimedia a,
    genetica.unidad b,
    genetica.ueb c
  WHERE a.iduni = b.iduni AND b.idueb = c.idueb;

for generate a layer on GeoServer but it needs an id field for linking to JSON and XML files. I created a table

CREATE TABLE genetica.gt_pk_metadata
(
  table_schema character varying(32) NOT NULL,
  table_name character varying(32) NOT NULL,
  pk_column character varying(32) NOT NULL,
  pk_column_idx integer,
  pk_policy character varying(32),
  pk_sequence character varying(64),
  CONSTRAINT pk_table PRIMARY KEY (table_schema, table_name, pk_column),
  CONSTRAINT unique_table UNIQUE (table_schema, table_name, pk_column)
)
WITH (
  OIDS=FALSE
);

and I filled with data about the view, but GeoServer does not work with that. Any suggestions?

I found this:
From my web search, the last message seems common among WFS-T users. So I clicked through several pages, tutorials and message boards (including this one) to find a solution. It seemed, that I needed a Primary Key Metadata table to get updatable views. So I added one with the name gt_pk_metadata_table to my database:

Then I referenced that table in my GeoServer data store (under "Primary key metadata table") and saved. Unfortunately, .....

Where do I create the Primary Key Metadata table?
How do I reference that table in my Geoserver data store?


texts: How can I convert my Shapely polygon to a numpy pixel array representing a boolean mask, with True for where the polygon is?


texts: Transferring values from field in table to field in attribute table of feature class using ModelBuilder?

texts: For example, I wanted to replace the values stored in the “enabled” field of P1 feature class with their corresponding values stored in “enabled-Sc” filed stored in J1. What tools might do this sort of thing effectively and efficiently within the capacity of model builder knowing that there is a common field between the stand-alone table (J1) and the attribute table (P1).


texts: Determining date or time Feature Class was created in File Geodatabase?

texts: I had a script running at the office over the last 5 days and tried to remote on unsuccessfully since the 3rd day. I got here and it appears there was a power outage. Is there a way to find out the time a feature class was created? That would very much help me figure out how long the script took to run.


texts: Getting incorrect coordinates for ortho photo using GDAL libraries?

texts:
I am very new to gis, and while most of the stuff I have is working well, this thing has me stumped.

I'm using GDAL libraries and the c# wrapper to access various gis functions. I'm processing ortho images to merge a few tiles, convert to WGS84 and then I chop the merged image into a bunch of small tiles.  For the most part, this is working well, but I've come across this image  which I can't get the right co-ordinate for.

PROJCS["NAD_1983_HARN_StatePlane_Virginia_North_FIPS_4501_Feet",
GEOGCS["NAD83(HARN)",
DATUM["NAD83_High_Accuracy_Reference_Network",
SPHEROID["GRS 1980",6378137,298.2572221010002,
AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","6152"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4152"]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",39.2],PARAMETER["standard_parallel_2",38.03333333333333],
PARAMETER["latitude_of_origin",37.66666666666666],
PARAMETER["central_meridian",-78.5],
PARAMETER["false_easting",37673535.76388889],
PARAMETER["false_northing",21527734.72222222],
UNIT["us_survey_feet",0.3048006096012192],
AUTHORITY["EPSG","2853"]]

I use the following code to extract co-ordinates (WGS84_STR is the WGS84 WKT):

string proj = ds.GetProjectionRef();
SpatialReference old_cs = new SpatialReference(proj);
SpatialReference new_cs = new SpatialReference(WGS84_STR);
trans = new CoordinateTransformation(old_cs, new_cs);

...then a method to return the coordinates:

...
ds.GetGeoTransform(adfGeoTransform);
dfGeoX = adfGeoTransform[0] + adfGeoTransform[1] * x + adfGeoTransform[2] * y;
dfGeoY = adfGeoTransform[3] + adfGeoTransform[4] * x + adfGeoTransform[5] * y;
trans.TransformPoint(adfGeoTransform,  dfGeoX, dfGeoY, 0);
...

where x and y are the image corner offsets.

I get lat of: -14.921860233519684, and lng of: -130.08927419763472 ...which is in the Pacific Ocean somewhere, not Northern Virginia. The tile is also rotated. I wonder if someone could point out where I've gone wrong?

Source of problem found
It turns out the GEOTIFF_CSV environmental variable was not set. In c#, the following code:
OSGeo.GDAL.Gdal.SetConfigOption("GEOTIFF_CSV", Path.GetDirectoryName(Application.ExecutablePath) + @"\gdal-data");
makes a very big difference. Now the geoTransform contains the actual lat/long.


texts: Using ApplyEdits in ArcGIS API for JavaScript?

texts:
In my ArcGIS Javascript 3.12 application I am drawing a freehand polygon, clicking a "Save" button and adding the polygon to my feature layer.  This works correctly except when the polygon crosses itself (Image 1).

When the polygon crosses itself and save is clicked, the success function runs and everything appears to be working.  When I open the attribute table in ArcMap to look at the polygon it has "empty" spatial values (Image 2).

Does anyone have any ideas?

I can provide code samples if needed.

Image 1

Image 2


texts: GeoNetwork multilingual entries on custom entry fields

texts:
We need the conformity pass to have three states, conformant(bool true), not conformant (bool false) and not evaluated (bool nil-reason). For this purpose we made the entry field for the conformity use a template that looks like this:

<template>

  <values>
    <key label="conformity_title"  xpath="" tooltip="gmd:pass">
      <helper name="gmd:title" context=""/>
    </key>

    <key label="conformity_date" xpath="" use="gn-date-picker">
      <directiveAttributes data-min-view="day"/>
    </key>

    <key label="explanation" xpath="" use="textarea"/>

    <key label="conformity_code" xpath=""/>

    <key label="conformity_codeSpace" xpath="">
      <helper name="gmd:codeSpace" context="gmd:CodeSpaceHelper"/>
    </key>

    <key label="pass" xpath="" use="gn-checkbox-with-nilreason">
      <directiveAttributes data-tag-name="gmd:pass"
     data-nilreason="eval#gmd:DQ_ConformanceResult/gmd:pass/@gco:nilReason"
     datalabels='{"true":"conformant","false":"notConformant","unknown":"notEvaluated"}'/>
    </key>
</values>

But if we use the template, we don't have multilingual entries for it.

Is there a way to either enable multilingual field or to make the stock entry field be a checkbox with nilReason?


texts: Merging multiple shapefiles and finding which shapefile gives the minimum value for each attribute in QGIS?

texts:
I have several different shapefiles that I want to combine into one based on the minimum value of each attribute. I found this solution. However, I also want to know which shapefile each minimum value came from.
For example, I have 5 shapefiles with values in each one of 10 attributes. For each attribute, I found the minimum value from one of the 5 shapefiles using the method in the solution above. But I also want to know which of the 5 shapefiles the minimum value came from for each attribute. Is this possible?


texts: ArcGIS 10.5 Python toolbox problem SelectLayerByAttributeERROR 000358: Invalid expression Failed to execute

texts:
I am trying to create a .pyt that find and select an object in an specific layer but I have a problem.
I'm working on the Desktop version 10.5 and actually the code is this one:

def getParameterInfo(self):
    """Define parameter definitions"""
    Cabina = arcpy.Parameter(
        displayName="Cabina",
        name="cab",
        datatype="string",
        parameterType="Required",
        direction="Input")

    params = [Cabina]
    return params

def isLicensed(self):
    """Set whether tool is licensed to execute."""
    return True

def updateParameters(self, parameters):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""
    return

def updateMessages(self, parameters):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""
    return

def execute(self, parameters, messages):
    """The source code of the tool."""
    mxd = arcpy.mapping.MapDocument("CURRENT")
    lyr = arcpy.mapping.Layer(r"ESERCIZIO\ELE_ESER_Cab")
    exp = "'" + str(Cabina) + "'"
    arcpy.SelectLayerByAttribute_management(in_layer_or_view= lyr, selection_type="NEW_SELECTION", where_clause= "COD_CABINA LIKE" + exp)
    arcpy.RefreshActiveView()

    return
    `

the problem is that continue to appear this error message: ERROR 000358: Invalid expression Failed to execute. (I'm working in an Oracle environment)
Actually I tried to put a value such as exp = 'A0001' and the script work well.
I read almost everything but seems like nothing works. I tried to find some documentation about how to use custom input and SQL but I did not find so much.


texts:
I have written a query that will create a table and perform a GROUP BY from an existing spatial table, but the table does have the geometry.  I am wondering where, and what, in this query I put to make the new table have spatial geometry (points).

Thanks...

SELECT st_id, name, COUNT(st_id)
    INTO climate.test_stats
FROM climate.lugdc_stations
GROUP BY st_id, name
ORDER BY count DESC;

ALTER TABLE climate.test_stats
ADD COLUMN pk_id serial PRIMARY KEY;


texts:
I've been looking for a while for Historical GIS datasets. In a lot of places (also in answers here at Stack Exchange), people are referring to a great dataset at ThinkQuest, which contains shapefiles for many years inbetween 2000 BC and 1994 for the whole world (I am interested in Europe, in particular).

But since ThinkQuest has been discontinued, the archive is only accessible through the ThinkQuest Library. Unfortunately, this archive does not contain the actual downloads of the datasets...

Does anybody know where I can find these datasets, or is there maybe someone who downloaded these in the past willing to share them here?


texts: What is the full syntax SQL reference for the `geometry` data type?

texts:
I have a database in which a field has the data type geometry(Geometry, 4326). In another similar database a field merely had the data type geometry.

Where can I find the full reference on this data type? Where can I find the list of possible values for the first argument to geometry()? I couldn't find them in the PostGIS manual.


texts:
I created kml file from ArcMap and imported into my Google maps, it works. Not sure if you can open this URL:
https://www.google.com/maps/d/edit?mid=13TuHmGuW9nnE-8TACTfGLdDDf6o&ll=38.90133009121122%2C-77.01476049999996&z=12

Now I want to display my map using Java API and I used this code
https://developers.google.com/maps/documentation/javascript/kml
tweaked it a little bit, inserted my API_Key at the end of the script but the problem is

<script>
      var map;
      var src = 'https://myGoogleMapURL

should be with the .kml extension if I am not mistaken.
My question is, how do I generate kml file/url from my map URL mentioned above.
I tried adding &output=kml  at the end of the Google map URL as some people suggested but it does not work.


texts:
Via the U.S. Census I have downloaded all of the census group shapefiles.  These are approximately 19-20GB across 56 shapefiles.

I want to merge these to do some testing.  Clearly, shape file format is not going to work, so I am using spatialite.  I have read that spatialite journalling causes performance issues when trying to load lots of data.  Two questions:

Can I turn off journaling and just bulk load the data.  I'll build stats and indices later.
Can I merge multiple shapefiles into one table within spatialite?  I know I could use OGR and write a script to pump the data into a table one shape file at a time, but is this possible direct via spatial lite.

Disclaimer:

I know that this will be an ~8,200,000 row table. I am not running this in production or planning on running many queries on this.  This is a research project.
I could also go with a full blown PostgreSQL install and have this running, but I want to dig a bit deeper into Spatialite.


texts: pgr_Dijkstra: iterate/loop through a table with multiple OD pairs one-by-one

texts:
I have a table with 2 columns: column A with the startvertex and column B with endvertex.

How can I do a OD pairs shortest path?

It is with a loop I suppose, but how can this loop be implemented into pgr_Dijkstra? It's impossible  do to all the OD pairs manually, so an iteration is necessary.

I have already a working query for pgr_Dijkstra for a single point-to-point, but the problem is to do it for all the OD-pairs one-by-one.


texts:
I have a desktop application for a client that needs to map a main subject location and around 5 "comparable" locations, within a range of 1-2 miles (often even less), from an aerial view, and zoom in as much as possible given the range of the locations. We have an interface to supply the locations, and are looking for software to simply generate and print the map.  No GUI needed, no fancy bells and whistles...

Can anyone recommend software to accomplish this?  I've looked at similar questions in many different forums, but the software recommended in the answers are usually overkill to the job I need to get done.

Thanks!

EDIT:
We've tried ArcGIS but it can't zoom as close as we'd like for aerial images, and doesn't print all areas correctly.


texts: Rectangle shapefiles with the same extents does not align perfectly after projection

texts:
I have a rectangle shapefile built with UTM19 coordinates using a python script, and then projected to NAD83_CSRS (EPSG:4617).

I have a second rectangle shapefile built with the same method (python with UTM19 coords), with the same corners coordinates but gridded with 200m squares. This layer is then also projected in NAD83_CSRS projection (EPSG:4617).

If I load the two layers in their original projection (UTM19), they align perfectly everywhere. But if I load the projected layer, They align perfectly in the corners but there is a shift when you look more in the middle area of the extents (see image). The shift is in the range of 5~8m.

I tried reprojecting the layers with ArcMap, QGIS and OGR and I always get the same result.

Anyone has an idea what that problem might be and how can I resolve it ?


texts: Keeping layer name when exporting to CAD?

texts:
I am trying to export selected features from layers  from SDE into a DXF using python, but it does not retain the layer name. But if I do it manually using the ArcToolbox tool, it does work.
Example below

Layers in MXD

DXF in AutoCAD  - which was created using the Export to CAD tool in ArcToolbox

DXF in AutoCAD - which was created using the python code below

The code that I run is below

mapDoc = r"C:\temp\DXF_Test10_v5.mxd"

mxd = arcpy.mapping.MapDocument(mapDoc)
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]

#Go through teh bookmarks layer
fields = "FULLNAME"
#location = "City of London"
location = "Golden Lane Estate"
expression = "\"FULLNAME\" = '" + location + "'"

for layer1 in arcpy.mapping.ListLayers(mxd):
    if layer1.name == "Bookmarks":
        with arcpy.da.SearchCursor(layer1, fields,expression) as cursor:
            for row in cursor:
                #select the record
                arcpy.SelectLayerByAttribute_management(layer1, "NEW_SELECTION", expression)
                #do select by location, select all the layers apart from bookmarks
                inFeatures = ''
                listOfLayers = []
                #loop through all the layers in teh map
                for layer in arcpy.mapping.ListLayers(mxd):
                    if layer.name <> "Bookmarks":
                        if layer.isGroupLayer == False:
                            print layer.longName
                            if layer.visible == True:
                                listOfLayers.append(layer)
                                #select the features that intersect
                                arcpy.SelectLayerByLocation_management(layer, "INTERSECT", layer1, "", "ADD_TO_SELECTION")
                                count = int(arcpy.GetCount_management(layer).getOutput(0))
                                print count
            #print inFeatures
            #export to cad
            #arcpy.ExportCAD_conversion(listOfLayers, "DXF_R2010", r"H:\Corporate GIS and Web Mapping\Data_Repository\DXF_TEST" + os.sep + location + ".dxf", "Ignore_Filenames_in_Tables", "Overwrite_Existing_Files", "")
            arcpy.ExportCAD_conversion(listOfLayers, "DXF_R2010", r"C:\Temp\DXF_COL" + os.sep + location + ".dxf", "Ignore_Filenames_in_Tables", "Overwrite_Existing_Files", "")

del mxd

Does anyone know why when I run the script it does not honor the layer name?


texts:
I'm trying to add a marker to a MapBox GL JS map ...

Here you are my code

    <meta charset="utf-8">
    <title>Test MapBox</title>

    <link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.37.0/mapbox-gl.css' rel='stylesheet' />

The code works fine but no marker is shown....

Suggestions / examples?


texts: Cannot input data when using tools e.g erase, but adding data to map works?

texts:
None of the data is available to use when inputting for tools use, but it all is available when just adding it?

Basically alot of my data is loaded into arcgis. When adding new data to my current map, there is no problem. All my downloaded data is readily avalable, but when it comes to using the existing data loaded onto arcgis already, when I want to for e.g intersect two layers, I cannot input anything, because as soon as I click input features, the data just isn't there anymore. And this happens for all other tools as well.


texts: Use two lists as input to Spatial Join in ArcPy?

texts:
I am trying to determine a method for conducting a spatial join on features from two separate feature datasets based on the first three letters of each input feature (A01_road == A01_gps).

import arcpy
import os

arcpy.env.overwriteOutput = True
arcpy.env.workspace = r'E:\GPS\GPS.gdb'

rds = r'E:\Pathways_GPS\Pathways.gdb\road_lines'
wlk_pts = r'E:\Pathways_GPS\Pathways.gdb\Walking_points'
output = r'E:\Pathways_GPS\Pathways.gdb\road_join'

#Spatial Join Variables
join_operation = 'JOIN_ONE_TO_ONE'
join_type = 'KEEP_COMMON'
field_mapping = ''
match_option = 'WITHIN_A_DISTANCE'
search_radius = '15 meter'
name_suffix = "_400join"

arcpy.env.workspace = rds
road_list = arcpy.ListFeatureClasses()
for road in road_list:
    road_fc = road[:3]
    arcpy.env.workspace = wlk_pts
    walk_list = arcpy.ListFeatureClasses()
    for walk in walk_list:
        walk_fc=walk[:3]
        if road_fc == walk_fc:
            out_fc = (os.path.join(output,(join_fc [:3]+ name_suffix)))
            arcpy.SpatialJoin_analysis(target_fc,join_fc, out_fc, join_operation, join_type, '', match_option, search_radius,'')


texts: Looking for free program, similar to 3D Analyst, to visualize tree stands in 3D

texts:
I want to make a 3D stand visualization of trees in a specific area. I need something similar in functionality to ArcGIS 3D Analyst, but I don't have a license for that plugin anymore.

The only alternative I have found is the Stand Visualization System from the USDA Forest Service, which works quite well, but has no GIS functionality.

I want to visualize the height and clustering of geographically referenced trees. There have to be different models for the tree species (spruce, beech, ... ).

EDIT:
For all people searching for more tools
-> here is another posting on gis.stack


texts:
According to this answer I can read in a .lyr file via OGR as long as it was made with ArcGIS version 10.

I have such a file (from the ArcGIS Maps & Data page) but I can't quite make R's readOGR function work with it.

Directory contents after decompressing the .lpk file:

[1] "ESRI_Zip_Codes_Dec2011.lpk" "iteminfo.xml"               "item.pkinfo"                "thumbnail.png"               [5] "USA_Zip_Codes.lyr"          "zip_poly.prj"               "zip_poly.sdc"               "zip_poly.sdc.prj"            [9] "zip_poly.sdc.xml"           "zip_poly.sdi"

And then trying to open the various components:

    ogrInfo( dsn=mypath, layer="zip_poly" )

  Error in ogrInfo(dsn = file.path(shapedir, "zipcodes", "ESRI_Zipcodes"),  :
   Cannot open file

    ogrInfo( dsn=mypath, layer="USA_Zip_Codes" )

  Error in ogrInfo(dsn = file.path(shapedir, "zipcodes", "ESRI_Zipcodes"),  :
   Cannot open file

How do I open these things?  I would normally take the easy way out and convert to Shapefile it via ArcGIS, but my annual license expired and my new one hasn't arrived yet.


texts: WARNING 000632: datum conflict when creating buffer in ArcGIS Desktop?

texts:
A buffer I created successfully executed, but I received the following message in ArcGIS Desktop 10:

  WARNING 000632: Datum conflict between map and output.

How can I fix this?


texts: Reprojection of nc file with gdal.Warp gives a NoneType object

texts:
In the following code I try to make a file which should have the same grid and bounds as the mask.

mask=gdal.Open(path1,gdal.GA_ReadOnly)
print(mask.GetGeoTransform())
src=gdal.Open(path2,gdal.GA_ReadOnly)
dst=gdal.GetDriverByName('NetCDF').Create(pathdst,xsize=mask.RasterXSize,ysize=mask.RasterYSize)
dst.SetGeoTransform(mask.GetGeoTransform())
print (dst.GetGeoTransform())
dst.SetProjection(mask.GetProjection())
bnds=(mask.GetGeoTransform()[0],mask.GetGeoTransform()[3]-mask.RasterYSize,mask.GetGeoTransform()[0]+mask.RasterXSize,mask.GetGeoTransform()[3])
optt=gdal.WarpOptions(format='NetCDF',xRes=mask.RasterXSize,yRes=mask.RasterYSize,outputBounds=bnds,width=mask.GetGeoTransform()[1],dstSRS=mask.GetProjection())
dst=gdal.Warp(dst,src,options=optt)
dst=None
src=None
mask=None

but it throws an error

  TypeError: in method 'wrapper_GDALWarpDestDS', argument 4 of type >'GDALWarpAppOptions *'`

Code

mask=gdal.Open(path1,gdal.GA_ReadOnly)
src=gdal.Open(path2,gdal.GA_ReadOnly)
bnds=(mask.GetGeoTransform()[0],mask.GetGeoTransform()[3]-mask.RasterYSize,mask.GetGeoTransform()[0]+mask.RasterXSize,mask.GetGeoTransform()[3])
optt=gdal.WarpOptions(format='netCDF',xRes=mask.RasterXSize,yRes=mask.RasterYSize,outputBounds=bnds,width=mask.GetGeoTransform()[1])
dst=gdal.Warp(pathdst,src,options=optt)
print(dst.RasterCount)
dst=None
src=None
mask=None

generates dst=NoneType and no file.

The code

dst=gdal.GetDriverByName('netCDF').Create(ppath+'test1.nc',mask.RasterXSize,mask.RasterYSize,1,gdalconst.GDT_Float32)
dst.SetGeoTransform(mask.GetGeoTransform())
dst.SetProjection(mask.GetProjection())
ds=gdal.ReprojectImage(src_prec, dst, prec_info[0],mask_info[0], gdalconst.GRA_Bilinear)
del ds,dst

works properly, but I can't do the same thing for multiband raster.

How can I reproject multiband src to mask domain with GDAL?

Today I tried to do the same thing in QGIS and python with same files. I took 1 band from src and saved it as nc. In QGis it worked perfectly:

gdalwarp -ot Float32 -t_srs EPSG:4326 -r bilinear -of GTiff -tr 1 1 -te 19.5 " 40.0" " 70.5" " 75.0" -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 -wo OPTIMIZE_SIZE=TRUE dst.nc

I copied all parameters to python and got NoneType result.

mask=gdal.Open(path1,gdal.GA_ReadOnly)#.GetRasterBand(1)
print(mask.RasterCount)
src=gdal.Open(os.getcwd()+'/test1/mask.nc',gdal.GA_ReadOnly)
bnds=(mask.GetGeoTransform()[0],mask.GetGeoTransform()[3]-mask.RasterYSize,mask.GetGeoTransform()[0]+mask.RasterXSize,mask.GetGeoTransform()[3])
optt=gdal.WarpOptions(format='netCDF',resampleAlg='bilinear',warpOptions=('OPTIMIZE_SIZE=TRUE'),creationOptions=('COMPRESS=DEFLATE','PREDICTOR=1','ZLEVEL=6'),outputType=gdalconst.GDT_Float32,srcSRS='EPSG:4326',dstSRS='EPSG:4326',xRes=mask.RasterXSize,yRes=mask.RasterYSize,outputBounds=bnds,width=mask.GetGeoTransform()[1])
dst=gdal.Warp(pathdst,src,options=optt)

gdal version -- 3.0.4


texts: I am new to Google Earth Engine and was messing around with the settings but did not find an option to switch the code editor to dark mode. I did a Google search and found a switch to dark mode for Google devtools but I don't believe that it applies to Earth Engine code editor. I also looked for a dark mode stylus on github and couldn't find one for Earth Engine. Does anyone know a simple way to switch to dark mode, or has a stylus I could use?


texts: Transfering styles from QGIS to ArcGIS Desktop?

texts:
We have a QGIS project that has symbols for a lot of different lines. Data comes from a geodatabase that is accessible to both QGIS and ArcGIS.

We now want to use the symbology of the QGIS project to style the features in the ArcGIS map document.

My first approach was to use SLDs. The QGIS style can be exported to SLD easily, but how can you style your features according to a SLD using ArcGIS Desktop?

Setting up a simple WMS does not suffice, both QGIS and ArcGIS Desktop need editing capabilities on the features.

If one could provide a SLD with a WFS (in a Esri-readable way) that would also be a possibility. Have found nothing for that, either.


texts: Importing csv file in QGIS including names and description

texts:
I'm just getting to learn QGIS. I want to make a map (Google maps) where I can map all our customers from our database. Including names and even better: description, address, telephone, etc. All this data is available in our own database.

So far I managed to import a csv list using Mmqgis.
Then added another layer with Openlayers
All works fine.

Also tried exporting the .shp as .kml, uploaded it to Dropbox and opened the link in a newly created personal Google Map. Works fine.

But the result is a map with locations. In which step can I include company names, not having to do this by hand...? I tried importing csv files via mmqgis with additional columns for names, etc. but I can't see where this information is stored in the QGIS project/I don't know if it is stored.

Any help is welcome!

Edit:
I found out that I can just import an xlsx file and select the correct columns in the google map creator...
Edit2: Now I see only Google accounts can edit the map, it would make more sense to use OSM or whatsoever I guess. Any thoughts on this? I'm diving into OSM now, but importing xlsx doesnt seem to be an option as easy as in Google Maps. So first question remains: how to include names/description in QGIS/OSM
Edit3: I'm looking for a solution csv->map including names/descriptions, preferably not Google Maps.


texts: How to upload vector data to OpenStreetMap ( quickly )?

texts:
I have 3 shapefiles containing:

8791 Linestrings
6704 Linestrings
584 Polygons

I want to upload this data to OSM ( share it ) I took a look at the guidelines of uploading data to OSM. And I want to know how should I structure my data in order to be conform to OSM norms.

So what I'm asking for the structure that I need to get my shapefiles in to be accepted ( attributes, types, ...).


texts: Writing and testing python code using existing OSGeo4W Python3 install?

texts:
Looking for a development environment, like idle, for Python 3.6.1 that I can install and use with my existing python3 installation from the OSGeo4W. I am trying to learn python for use with QGIS and don't want multiple python installs of various versions confusing me and my work environments.

How should I set myself up for best results?


texts: Merging csvs to shapefiles and plotting county maps using time series loop in R

texts:
I have eleven county scale shapefiles named "df_1_imid_soy_2004.shp" for all years from 2004 to 2014 containing rows for every county in three states (Dakotas and Minnesota)

I also have eleven csv files containing a "st_cnty" unique key that is used to merge to its corresponding shapefile. There is also a column called "g_per_acre_soy" that is the continuous variable I want to plot with a yellow to red color map to illustrate the quantity.

library(ggmap)
library(rgdal)
library(rgeos)
library(maptools)
library(dplyr)
library(tidyr)
require(sp)
require(RColorBrewer)

par(mfrow=c(4,3),mar=c(0.6,0.6,0.6,0.6), oma=c(1,1,1,1))

yearList = c(2004:2014)
for (year in yearList){
  #shape is original shapefile and table is csv to join
  shape = readOGR(".",print(paste0("df_1_imid_corn_",year)))
  table = read.csv(print(paste0("df_1_imid_corn_",year,".csv")))
  shape@data$order <- 1:nrow(shape@data)
  shape@data <- merge(shape@data, table, by.x = "st_cnty",by.y="st_cnty",all.x=T)

  #new shapefile to remove nas because not all rows join
  shape2 <- shape[!is.na(shape@data$g_per_acre_corn),]
  shape2@data$order <- 1:nrow(shape2@data)

  #the variable being plotted has a large variance between counties.
  #here I attempt using breakpoints and cut and assign each row a color
  breakpoints <- seq(0,1,by=0.01)
  myPaletteBlue = brewer.pal(1,'RdYlGn')
  createColors  = colorRampPalette(myPaletteBlue)
  myColorsBlue  <- createColors(length(breakpoints)-1)
  numberColors = length(unique(myColorsBlue))
  shape2@data$g_per_acre_corn_cat = cut(shape2@data$g_per_acre_corn,breakpoints)
  shape2@data = shape2@data[order(shape2@data$g_per_acre_corn),]
  colorOfCats  = data.frame(g_corn_cat=levels(shape2@data$g_per_acre_corn_cat),colors=unique(myColorsBlue))
  shape2@data = merge(shape2@data,colorOfCats)

  #then make the plots
  plot(shape2,col=as.character(shape2@data$colors))
  plot(shape, add=T)
  title(main=year,add=T)
}

The merging and plotting works, but I'm not sure how best to plot the column with high variance and significantly increases from year to year. For example, the dependent variable in 2004 ranges from (0.1-250,n=110,mean=2) but in 2014 ranges from (0.9-500,n=130, mean=50).

So I would like to plot all years on the same figure with ONE INDIVIDUAL SCALE to visualize annual changes.

1) What is the best way to keep the same scale for all years, but also have enough variation with colors to visualize inter-county changes within the same year?

Below, red=data, white=NA (not joined)


texts: Geonode layers not visible on map creation

texts:
On a geonode 2.4 installations I have hundreds of layers. When creating a new map (as logged in admin) I can see none of them in the available layer tree.  It just stays empty. geonode.log does not show any error.

I recognized that the Geoserver which is running on same machine is not in the list.

Thus, I tried to add it like http://mygeonode.org/wms/ unfortunately this does not work.

This is how my server definition on local settings looks like:

    GEOSERVER_URL = SITEURL + 'geoserver/'

# OGC (WMS/WFS/WCS) Server Settings
OGC_SERVER = {
    'default' : {
        'BACKEND' : 'geonode.geoserver',
        'LOCATION' : 'http://localhost:8080/geoserver/',
        'PUBLIC_LOCATION' : GEOSERVER_URL,
        'USER' : 'admin',
        'PASSWORD' : 'geoserver',#'password',
        'MAPFISH_PRINT_ENABLED' : True,
        'PRINT_NG_ENABLED' : True,
        'GEONODE_SECURITY_ENABLED' :True,
        'GEOGIG_ENABLED' : False,
        'WMST_ENABLED' : False,
        'BACKEND_WRITE_ENABLED': True,
        'WPS_ENABLED' : True,
        'LOG_FILE':'/usr/share/geoserver/data/logs/geoserver.log',
        # Set to name of database in DATABASES dictionary to enable
        'DATASTORE': 'datastore',
    }

What´s wrong with my local geoserver configuration?


texts:
My question is a follow-up of this:
plotting US cities on a map with matplotlib and basemap.
I would like to use some Python library to plot on a map city locations, starting from a list like this:
    L = ['Paris', 'Lyon', 'Marseille', '...']

In the linked GIS question, this was easily achieved through Basemap, which has however been deprecated in favour of Cartopy and does not currently work on Pyhton 3. Is there a simple way to achieve the same result with Cartopy or another library?
This Reddit thread may be helpful for advanced users but being a newbie I was looking for something more intuitive.


texts: Publish a Query Layer as a Feature Access Service - ArcServer

texts:
I am connecting to a SQL Server 2016 database, and adding data into ArcMap 10.3.1 as a Query Layer.

I want to publish this data as a Feature Access service to ArcServer 10.3.1. I have no problems publishing this data as a mapping service. However, when I want to enable Feature Access, I get 00223: Query layer has virtual column error.

I have read the documentation on this error, and the only solution is to "Change the query definition of your layer so it no longer contains a virtual column or columns".

I have tried several different ways to alter the sql query (very simple queries), and all still produce this same error.

What is a virtual column? How can I publish a Query Layer as a Feature Access service?


texts:
I have been using ArcGIS and after sometimes suddenly it has stopped working.
When I tried to open ArcMap, it brings an error message:

ArcGIS has stopped working. Please wait while gathering information about the error

then  nothing will appear afterwards.
I am using 10.7 version.
I have tried to reinstall, repair, but failed to work.
Has experienced the same error?


texts: I might be thinking about this all wrong, but right now I'm using Mapnik's OGC Server, and using that WMS with OpenLayers. How would I group multiple layers (say country, water, state borders) into one layer, so that all of those combined show up as a Base Layer in OpenLayers?  I'm trying to do this through the Python API.


texts: Some labels not being generated on map contours

texts:
I'm trying to produce a map using the latest QGIS incarnation: 2.18.11
After importing the UK OS contour data for one 100km grid square, I crop it using two contiguous polygons (1km squares). Two of the contour lines become cut into two separate parts. When applying labelling to the contour lines layer, I find that only one part of the two part lines displays the labels (50m & 60m in the attached). The two parts of the two lines have the same ID. Any ideas anyone?


texts: Is it possible to have auto-updating area fields?

texts:
I'm using QGIS 2.4.0 and I've noticed some problems with the derived area and the $area field.
Firstly, I have a polygon which is a filled ring, inside another polygon. The derived area, and therefore my $area field, is saying 600sqm, but the actual measurement made with the area measure tool, is 300sqm. Why the discrepancy?

Also, why am I having to re-write the area field every time I make any changes to the vector layer. It's a real pain as I'm having to present areas in Hecatre so the formula involve both conversion and rounding. To have to re-write it each time I make any adjustment is really slowing me down. Any plugins that might have an auto-updating area field?


texts: How to extrude a single polygon shape to a interpolated surface elevation in arcscene?

texts: I have a new problem to solve which is: how to extrude a single polygon shape to the curvature of a defined surface but not only to a specific one value height? I have several interpolated surfaces and I want to extrude a rectangular polygon shape inclined to those surfaces so that I can make a block diagram of several geological layers.:-)


texts:
I would like to create a raster mosaic of the national USGS NED 30m DEM (~50 GB). I'd like to leverage parallel processing using python multiprocessing module where possible. I'm using ArcGIS 10.0. My first approach was divide and conquer - mosaic 2-10 at a time and then mosaic the result.

Q) What is the best method for mosaic parallization? Examples?

EDIT: the original two-part question was split into two questions. The other is here.


texts: Exporting images from a list to drive fails without error in Google Earth Engine

texts:
I am mapping over a function that returns images, more specifically changemaps made with the LandTrendr API, which are ee.Images (according to the LandTrendr for GEE homepage).
I end up with a list (see print output in code below) of images and want to use a for-loop to write them to my Google drive. (I am aware that for-loops are deprecated in GEE.)
However, no task is created and no error is thrown, leaving me puzzled to what is wrong/missing.
Here is a reproducible example. If you have not used LandTrendr before and want to run this example you have to visit this site first to add the API to your account.
var years = ee.List.sequence(1985, 1995, 5);
var aoi = ee.Geometry.Point(13.64, 55.9552).buffer(500);

var imageList = years.map(function(n) {

// define landtrendr parameters
var startYear = 1985;
var endYear =  1995;
var startDay = '06-20';
var endDay = '09-10';
var index = 'NDVI';
var maskThese = ['cloud', 'shadow', 'snow', 'water'];

var runParams = {
  maxSegments:            6,
  spikeThreshold:         0.9,
  vertexCountOvershoot:   3,
  preventOneYearRecovery: false,
  recoveryThreshold:      0.9,
  pvalThreshold:          0.05,
  bestModelProportion:    0.75,
  minObservationsNeeded:  6
};

var strt = ee.Number(n)
var nd = ee.Number(n).add(5)

// define change parameters
var changeParams = {
  delta:  'loss',
  sort:   'greatest',
  year:   {checked:true, start: strt, end: nd},
  mag:    {checked:false, value:100,  operator:'>'},
  dur:    {checked:false, value:5,    operator:'<'},
  preval: {checked:false, value:300,  operator:'>'},
  mmu:    {checked:false, value:11},
};

// load the LandTrendr.js module
var ltgee = require('users/emaprlab/public:Modules/LandTrendr.js');

// add index to changeParams object
changeParams.index = index;

// run landtrendr
var ftvList = ['NDVI'];
var lt = ltgee.runLT(startYear, endYear, startDay, endDay, aoi, index, ftvList, runParams, maskThese);

// get the change map layers
var changeImg = ltgee.getChangeMap(lt, changeParams);

return changeImg

});

print(imageList) //Output below

//List (3 elements)
  //0: Image (6 bands)
    //type: Image
    //bands: List (6 elements)
      //0: "yod", signed int16, EPSG:4326
      //1: "mag", double, EPSG:4326
      //2: "dur", double, EPSG:4326
      //3: "preval", double, EPSG:4326
      //4: "rate", double, EPSG:4326
      //5: "dsnr", double, EPSG:4326
  //1: Image (6 bands)
  //2: Image (6 bands)

for (var i = 0; i<imageList.length; i++) {

  // get changemap from list
  var im = imageList.get(i)
  // Export band named 'mag' as GeoTiff
  Export.image.toDrive({
  image: im.select(['mag']),
  description: 'changemap_'+i,
  scale: 30,
  region: aoi,
  fileFormat: 'GeoTIFF',
  maxPixels: 3991919404

});

}


texts: Creating Raster with values based on function output

texts:
I have a function that calculates a price based on the input slope and distance. I want to write the price to a raster as the rastervalue. How do I do that? OpenSource and ArcMap solutions would work.

slopeRaster = "slope.tif"
emptyRaster = "emptyraster.tif" # How do I create an empty raster?
road = "road.shp"

for cell in emptyraster:
    # get slope from sloperaster at cell location
    ...
    slope = ...

    # get distance to nearest road from center of cell
    ...
    distance = ...

    # calculate price for cell
    price = pricefunc(slope, distance)

    # write price to cell as value  # How do I write a value to a raster


texts: How to select features inside polygon instead of inside polygon's bbox(extent) in openlayers 3

texts:
ol3 supplies the method forEachFeatureIntersectingExtent which only good for boxes.

I have user drawn polygons (or multipolygons) for selecting features from a vector layer.

I'm relatively new to openlayers and chose to start developing a new project with ol3.

I know openlayers 2 had intersect functions for geometries not just for extents.

How do I achieve this with ol3?


texts: How to specify raster name and how to designate the saved name of a created raster

texts:
import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")
env.workspace = r"G:\GHY_568\project"

mystery_creek = r"G:\GHY_568\project\Mystery_Creek.shp"
riofajardo_stations = r"G:\GHY_568\project\riofajardo_stations.shp"
gauge_station = r"G:\GHY_568\project\whiteoak_gauge.shp"
dems = arcpy.ListDatasets("*_filled", "Raster")

#print dems
#bethel_filled = r"G:\GHY_568\bethel_filled"
#arcpy.CheckOutExtension("Spatial")

for dem in dems:
    arcpy.env.extent = dem
    arcpy.env.cellSize = dem
    out_flow_direc = FlowDirection(dem)
    out_flow_acc = FlowAccumulation(dem)
##    if out_flow_acc > 200:
##        dem = 1
##    else:
##        dem = 0
    out_flow_length = FlowLength(out_flow_direc)
    out_cost_path = CostPath(gauge_station, out_flow_direc, out_flow_direc)
    if dem = "bethel_filled":
        out_watershed = Watershed(out_flow_direc, mystery_creek)
    else:
        out_watershed = Watershed(out_flow_direc, riofajardo_stations)

I need to specify which dem to use in the if statement and wasn't sure how to handle that, with "bethel_filled" in quotes or the path name or whath. Also the if statement at the end is in order to choose which dem to preform the watershed on because they have two different pourpoints. How do I specify which out_flow_direction to use in the if else statement. Mainly, how do I designate the saved name of each out_flow_direc.


texts: How to convert GeoSeries polygons to lists of coordinates?

texts: I have a GeoDataFrame consisting polygons. How can I convert each polygon into a list of coordinates? Is there any attribute, method or function in GeoPandas to do this?


texts: Download Dataset From ArcGIS REST Service

texts:
Looking at this information page for a dataset hosted on esri.com:
http://fema-services2.esri.com/arcgis/rest/services/2012_Sandy/ImageCat_NLT/MapServer/layers

Trying to figure out how I can get access to the raw data (either raw lat/lng coordinates or SHP).

It looks like the data is available, I just can't figure out how to pull it from this interface.


texts:
I have some mining data that I want to filter by year and by bench. The data is something like this:
bench   year
1      2020
1      2021
1      2022
2      2020
2      2021
2      2022
4      2020
4      2020
4      2021
4      2021
4      2022
3      2020
3      2021
3      2022
3      2023

This is a polygon layer with squares representing each bench. Each bench has a different year when they will be mined. some of these polygons have the same location (same X,Y but different Z value), different bench number, but will be mined at a different year.
I want to create several maps that can filter the data by year; e.g. These are all the benches that are going to be mined in 2020. Then I want to do it by bench; e.g bench 4 will be mined in 2020, 2021 and 2022.
I was wondering if I can do this with Atlas?


texts: OpenLayers Plugin QGIS Essen 2.14.3 defect?

texts:
I'm using QGIS Essen 2.14.3 for analysing gpx-routes.

After installing the latest Version, it is not possible to install "OpenLayers-Plugin". Python says: "module" object has no attribute "wkbCurve".

Can this error be avoided?

I need OMS, because I want to show only buildings, ways in the background of the gpx-routes.


texts: Editing curved annotation - unable to add vertex to baseline sketch, move bezier curve handles independently

texts:
It appears as though I am unable to add a vertex to a curved annotation during an editing session in ArcMap. When the annotation is right-clicked, the box to add or delete a vertex appears, but the options are greyed out.

While also on the curved annotation editing question - is it possible to move the bezier curve handles independently of one another? When I move one, the opposite handle moves as well and changes the baseline on both sides of the vertex.

It would be quicker to just re-do the smaller annotations, but sometimes my job involves editing long, meandering text strings and it is nice to avoid awkward gaps that result from certain curves.

edit - It appears that it is only the annotations in this specific map document, as the annotations can be edited in any other .mxds which means this problem will probably never get get figured out.


texts: Converting GeoPDF to vector format like CAD or shapefile?

texts:
I want to convert a GeoPDF that has vectors entity in it to a vector format and preserve the real coordinates of these vectors. The GeoPDF has been exported from ArcMap and has layers in it.

I have try GDAL, but it rasterize the PDF before exporting it. Which is not good for me. Also I don't want a solution that would convert this image back to vectors, too complicated for the data I have.

Is someone aware of any tool that can do that? The resulting format can be a CAD file or a shapefile or anything that is a vector format.


texts: Join attributes of point and line layer by location (closest distance) in QGIS

texts:
My problem is the following: i have a line shp (target layer) and a point shp (join layer). the second one (point layer) holds an attribute which i want to join with a line layer based on closest distance to the line. Points and lines are not intersecting

I can do this with ArcGIS and spatial join(match option: closest), but i don't know how can I achieve this in QGIS.


texts: Changin features' pointRadius based on the current zoom level

texts:
I need to have my features change their pointRadius based on the current Zoom level. My questions are:

Would it be possible to use Rule Based Styling techniques to overcome the problem?
If yes, it will work on map zoom changed as well?
If no, what are the alternative solutions and what you offer instead? (except redrawing features)


texts:
I am getting an error when using map.addLayer or marker.addTo(map) it gives

  t.onAdd is not a function

I believe either I am missing any library or what.
please help me out to get rid of this problem.
latLongShipArray contains lat/long values

var array = [];  // create an global array

for (var i = 0; i < latLongShipArray.length - 1; i++) {  //latLongShipArray consists lat/long values

    if (shipCoordinateArray[i]["Type"][0] == 1) { // condition to add marker ,shipCoordinateArray have some conditions
        marker = new L.marker(latLongShipArray[i], { //add lat/long in new marker
            icon: shipIConBigger
        });
        //marker.addTo(map);
        //map.addLayer(marker);
        marker.on('click', onMarkerPointClick);
    }
    else {   //condition to add marker
        marker = new L.marker(latLongShipArray[i], { // add lat/long in new marker
            icon: shipIConSmaller
        });
        //marker.addTo(map);
        //map.addLayer(marker);
        marker.on('click', onMarkerPointClick);  //click event
    }
    array.push(marker);   //push marker to array
}
var shipCoordinateMarker = L.layerGroup([array]);  // add a layergroup
//shipCoordinateMarker.addTo(map);
map.addLayer(shipCoordinateMarker);  // add marker array to map


texts: I am trying to import .mpk file in ArcMap 10.3.1. But, when I try to open using folder connection, the file with .mpk extension doesn't show up.. It only shows .tif file but not mpk file.


texts:
I need to do the following task: to retrieve the road type (highway, bridge, track, pedestrian and others) by coordinates, query a remote service (via HTTP).

I'm looking at the different major solutions available on the web (OpenStreetMap, Google Maps, Bing Maps, HERE Maps, Apple Maps), trying to compare them to better identify their limitations, but I need help completing this list.

At the moment, I just scouted OpenStreetMap service:
using OverpassAPI or a 2-queries method: using Nominatim + XAPI

What about the other services?

Google Maps: ???
Bing Maps: ???
HERE Maps: ???
Apple Maps: ???


texts: Setting coordinates of map centre from featureCollection geometry?

texts:
I would like to set the coordinates of the map centre using the coordinates from a geometry.

// randomly select 100 points
points = points.randomColumn('random');
var P = ee.Number(100).divide(points.size());
points = points.filter(ee.Filter.lt('random', P));
print (points)

// plot them
Map.addLayer(points);
var centre = points.geometry().centroid().coordinates();
print (ee.Number(centre.get(0)), ee.Number(centre.get(1)));
Map.setCenter(ee.Number(centre.get(0)), ee.Number(centre.get(1)), 9);
// Map.setCenter(-0.06657026751372405, 51.47080255078299, 9);

I am getting the error

  Provided center object has invalid values for lat, lon, or zoom.
  Expected numeric values.

Link to example


texts: Not able to save edited layer connected via PostGIS in QGIS

texts:
I have project and I was doing editing layers in QGIS. These layers are connected in QGIS via PostGIS. One layer of the four layers was saved correctly but other layers gave me error message.

Error message looks like this:

  ERROR: 15 geometries not changed. ERROR:  LineString must have at least
  two points CONTEXT:  SQL function "st_geomfromwkb" statement 1


texts: From SQL database to pandas dataframe to numpy array to ESRI Table

texts:
End goal: Join data pulled from a SQL database to three different feature classes in a ESRI geodatabase.

Efforts:

Using PYODBC, I've connected to the database and dumped the data into a Pandas Dataframe:

import pandas as pd
import pyodbc
import arcpy
import numpy as np

cnctStr='Driver={SQL Server Native Client 11.0};Server=XXXXX;Database=XXXXX;Trusted_Connection=yes;'
conn = pyodbc.connect(cnctStr)

SQL_Query = pd.read_sql_query(
'''SELECT        FILE_NUMBER, TYPE_CODE, STATUS_CODE, REVIEWER_ID, PREVIOUS_REVIEWER, OWNERS, WATERBODY_NAME, PROJECT_CATEGORY, PRELIMINARY_CATEGORY, APPLICANTS_REQUESTED_IMPACT, CALCULATED_AREA_OF_IMPACT, IMPACT_WETLAND_TYPE_CODE, CITY_NAME, COUNTY, LAST_ACTION_DATE, DATE_LAST_INFO_LETTER_SENT, IS_MORE_INFO_LETTER_SENT, DATE_MORE_INFO_LETTER_RETURNED, DATE_RECEIVED
FROM            WETLANDS_PERMIT_V
WHERE        (STATUS_CODE IN ('REVIEW EXT', 'PENDING', 'EA_ISSUED', 'MI RECEIVE', 'MI REQUEST')) AND (DATE_RECEIVED >= '01/01/2014')''', conn)

df = pd.DataFrame(SQL_Query, columns=['FILE_NUMBER','TYPE_CODE','STATUS_CODE','REVIEWER_ID','PREVIOUS_REVIEWER','OWNERS',
                                      'WATERBODY_NAME','PROJECT_CATEGORY','PRELIMINARY_CATEGORY','APPLICANTS_REQUESTED_IMPACT',
                                      'CALCULATED_AREA_OF_IMPACT','IMPACT_WETLAND_TYPE_CODE','CITY_NAME','COUNTY','LAST_ACTION_DATE',
                                      'DATE_LAST_INFO_LETTER_SENT','IS_MORE_INFO_LETTER_SENT','DATE_MORE_INFO_LETTER_RETURNED','DATE_RECEIVED'])

Which works beautifully.  But, to get the dataframe into ArcGIS, I have to get it into a numpy array format because the appropriate tool arcpy.da.NumPyArrayToTable expects a numpy array.

Using the new pd.to_numpy() leads to this error:

nparray = df.to_numpy()
arcpy.da.NumPyArrayToTable(nparray, r'C:\Temp\All_Wetland_Permits_DES_EA_Dash.gdb\Wetlands_Permit_V_Output')

**TypeError: narray.fields require**

Fix attempt 1:

Found advice here: https://community.esri.com/thread/226531-arcpydanumpyarraytotable-returns-runtimeerror-create-table
to try and use this code:

x = df.reset_index()
z = np.rec.fromrecords(x, names=x.columns.tolist())
arcpy.da.NumPyArrayToTable(z, out_table)

but that leads to this error:

TypeError: 'in_array' is not numpy array


texts: How to convert Raster image to Shapefile WITHOUT flipping it upside down

texts:
I am using QGIS v 3.2.1 to convert Raster image to Shapefile

Unfortunately the resulting shapefile is upside down.

How to keep resulting shapefile right side up

Here is animated gif (link on GitHub)


texts: Connect coordinates only according to specifications in an Excel list

texts:
I have an Excel list with coordinates of shops and their names. I want to connect certain coordinates with the help of a second list.

The second list contains in the first row and column the names of the shops. If there is a visual relationship between the shops, the list contains a "1" otherwise the field is empty.

Is there a possibility in QGIS to connect only those coordinates, where a "1" is entered in the second list?