texts: Dealing with NAN values in my training dataset

texts:
I am new to remote sensing. I am working on a project where I need to perform land use land cover classification. For such a task I want to use machine learning (python).
My features contain 'BANDS' (6 bands), 'NDVI', 'NDWI', and 'NORM' extracted from the sentinel level 1C images. I did cloud masking and then interpolation (to fill the gaps) and finally, I sampled the data and split them into training and testing datasets. I have also used a mask to eliminate samples where I have '0' which essentially means 'no-data':

# remove points with no reference from training (so we dont train to recognize "no data")
mask_train     = labels_train == 0
features_train = features_train[~mask_train]
labels_train   = labels_train[~mask_train]

# remove points with no reference from test (so we dont validate on "no data", which doesn't make sense)
mask_test     = labels_test == 0
features_test = features_test[~mask_test]
labels_test   = labels_test[~mask_test]

However, to be extra sure, I checked for the NAN values using:
 np.isnan(train_test).any() (for both training and testing datasets) I found that my features had NAN values. Any idea how to solve this? and why do I have NAN values in the first place?

Edit:
This is what my dataset looks like ( a lot of NAN values):


texts: Which CRS would be the most appropriate to georeference a map of Egypt?

texts:
I'm trying to georeference a geological map of Egypt from 1981, but I don't know which CRS to chose, either the transformation type nor the resampling method...

And once the georeferencing of the map is made, is there any way to convert it to another CRS? For example, WGS 84 or another one more standard?

I'll let you the image of the map where the only information given is that it uses the Lambert Conformal Conic Projection.


texts: WFS how to consume with different coordinate system?

texts:
I have a WFS in EPSG:4167. My map view has EPSG:3857. How to make the WFS work on a map view with different coordinates?

I used this to define my coordinate system from wfs:

import { register } from 'ol/proj/proj4.js';
proj4.defs("EPSG:4167", "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs");

register(proj4);

Note it is a third party wfs, so I only can use this EPSG.

Any ideas?

My first thought is that I should get the geojson from WFS and then transform the entire geojson using proj4 lib and display it on map (I'm using openlayers). Not sure if this is the best approach.


texts: How to display Flow Direction arrows + Burn Stream Netwok in QGIS?

texts:
I have used the automated r.watershed tool in QGIS 2.18 (GRASS plugin) to conduct a watershed analysis (basin boundaries, flow direction, flow accumulation etc).

I'm stuck on the following operations:

1. Visualising flow direction arrows. I tried to convert the flow direction raster to a vector point file using GRASS>r.to-vect, but it keeps crashing (not responding/no error message). Please can anyone recommend a workflow for this (converting the file + visualising arrows)?

2. Burning stream network into the DEM. I created a polyline file from the flow accumulation raster, to join up 'missing' channels. After converting the polyline to a raster, the new raster has a smaller extent than the original DEM. Thus, the SAGA>Hydrology>'Burn stream network into DEM' tool won't operate with rasters of different extents.


texts: gdalwarp -wo SOURCE_EXTRA option, what does it do?

texts:
I see people using the -wo SOURCE_EXTRA=100 option in their gdalwarp command when reprojecting to WGS84. What does it do? In the documentation it says:

"-wo" Set a warp option. TheGDALWarpOptions::papszWarpOptionsdocs show all options.  Multiple-wooptions may be listed.

But I can't find this "TheGDALWarpOptions::papszWarpOptionsdocs" doc anywhere. Anyone able to explain or even better point me towards this document?


texts:
I have a POSTGIS database with geography points and I want to cluster the points using quad-tree. I have a map where I have to display the markers clustered or non clustered depending on the zoom level so I can't really use k-means algorithm because I can't fix k.

How can I implement quad-tree in my application?

Please tell me if I am wrong somewhere.

What I have tried
I have tried clustering in the front-end but it doesn't help to solve my problem.
With a large number of markers that doesn't work. I want to send the clusters instead of all the markers.


texts: Uploading Geocode Service Definition Continues To Fail

texts:
I am not sure why this is happening, I am able to stage a Geocode Service, but when I try to  upload it to the server it continues to fail. I am able to upload regular MXD Mapping Services.
Below is the error that I keep receiving, but there is nothing of note in the server logs.
My only guess is that the server is 10.2.2 and I am using ArcGIS Desktop 10.3, but I didn't think that mattered since 10.1


texts: Making thematic map for each column in table using ArcGIS Desktop?

texts:
I have a table with header like
LocationId, Jan, Feb, Mar,.....Dec

I would like to make a graduated symbol map for all the locationIDs for each month, and make an output as a jpg/pdf file. What is the easiest way to do this in ArcGIS?

I can work on it one by one, but what if I get 20 columns to work with? Do I need to do it using ArcObjects?


texts: Moving points in ArcPy given complex requirements (Looping through polygons and extracting data)

texts:
How do I loop through a set of polygons in arcpy, using those polygons to
extract other data such as points and other polygons? (which I can then run analysis on)


texts: Imported DMS points from CSV to a point layer. Do I need to convert DMS points to Decimal Degrees to use metric distance units for Buffer?

texts:
DMS - Degrees Minutes Seconds

Buffer tool only shows degrees as the distance unit (not meter or kilometer)

Do I need to convert DMS points to Decimal Degrees to use metric
distance units for Buffer?
Is there a plug-in that automatically converts DMS to DD in QGIS?

Or is


texts: Quickly toggling visibility of certain symbology on ArcGIS Pro?

texts:
I'm trying to migrate from QGIS to ArcGIS Pro (still on trial version) and I have really basic problem.

How to quick toggle visibility of object with certain symbology (unique values symbology)?

In QGIS there is option to on/of visibility for symbol of every unique attribute value on Layers Pane. However, in ArcGIS PRO there is no option like that.

Here is how it look on QGIS:

ArcGIS Pro Contents pane:


texts: Checking Points within polygon and assigning value using ArcPy

texts:
How can I write arcpy code to determine if the points of a point feature class are inside polygons of a polygon feature class and write the result (true/false) to the point FC's attribute field?

I need to use the cursor function, but I don't know where to use it.

# Import arcpy module
import arcpy

# Script arguments
Point = arcpy.GetParameterAsText(0)
if Point == '#' or not Point:
    Point = "Point" # provide a default value if unspecified

Polygon= arcpy.GetParameterAsText(1)
if Polygon== '#' or not Polygon:
    Polygon= "Polygon" # provide a default value if unspecified

# Local variables:
Point__2_ = Point
Point__3_ = Point
Point__5_ = Point__3_
Point__4_ = Point__3_

# Process: Add Field
arcpy.AddField_management(Point, "REST", "TEXT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

# Process: Calculate Field
arcpy.CalculateField_management(Point__3_, "REST", "\"False\"", "PYTHON", "")

# Process: Select Layer By Location
arcpy.SelectLayerByLocation_management(Point, "COMPLETELY_WITHIN", Polygon, "", "NEW_SELECTION", "NOT_INVERT")

# Process: Calculate Field (2)
arcpy.CalculateField_management(Point__2_, "REST", "\"True\"", "PYTHON", "")```


texts:
i want to ask about erase outside in MapBasic.
I have two table : sensus and voronoi
And I want to erase the overlap areas. How can I do it?

This is my syntax

Select overlap(sensus.Obj,voronoi.Obj)
from sensus, voronoi
where sensus.Obj intersects voronoi.Obj

but that's indicates error 'The object in argument1 cannot be a point'

This is my full source code

'Create Table "voronoi"
    (Block_no Char(15),Remark Char(10),Type_Palm Char(10),Ha Float)
    File "D:\voronoi.TAB"
    TYPE NATIVE Charset "WindowsLatin1"
Create Map For Voronoi CoordSys Earth Projection 1, 104

Set Table Voronoi FastEdit On Undo Off
Create Object As Voronoi from sensus Into Table voronoi
Update voronoi
    Set Ha = Area(obj, "sq m")
Commit Table Voronoi

Set Table Voronoi FastEdit On Undo Off
select * from Voronoi
    where Ha >= Val(inputnumber)
    into __TO__UPDATE NoSelect

Set target on
Select overlap(sensus.Obj,voronoi.Obj)
from sensus, voronoi
where sensus.Obj intersects voronoi.Obj

Commit Table Voronoi
Close Table __TO__UPDATE

Create Index On voronoi (Block_no)
Add Map Layer voronoi
select * from Voronoi
    where Ha >= Val(inputnumber)'

I want to eliminate the line outside the border. Please help me


texts: Position of tiles using MapProxy server is wrong

texts:
I configured and deployed Mapproxy server on Ubuntu-14.04 and trying to get TMS service from Mapproxy.
 When I view my map in a web browser by referencing the host server directly(MapProxy Demo), everything is seen appropriately. But, when the same request is processed via Leaflet webpage, the order of tiles  is messed up (the tile position is wrong.)

I thought something is wrong with projection so I tried to fix it using
Proj4Leaflet API (code included below.) What can be the appropriate projection and CRS to fix this in leaflet? Is there any other solution?

Following is my code in the .yaml file

services:
  demo:
  tms:
    use_grid_names: true
    # sets the tile origin to the north west corner, only works for
    # tileservice at /tiles. TMS at /tms/1.0.0/ will still use
    # south west as defined by the standard
    origin: 'nw'
layers:
  - name: tms
    title: TMS
    sources: [osm_cache_in]
caches:
  osm_cache_in:
    grids: [osm_grid]
    # cache will not be stored locally
    # disable_storage: true
    # a tile source you want to reproject
    sources: [osm_source]
sources:
  osm_source:
    type: tile
    grid: osm_grid
    url: http://a.tile.openstreetmap.org/%(z)s/%(x)s/%(y)s.png
grids:
  osm_grid:
    base: GLOBAL_MERCATOR
    srs: 'EPSG:3857'
    origin: nw
globals:
  # cache options
  cache:
    # where to store the cached images
    base_dir: './cache_data'
    # where to store lockfiles for concurrent_requests
    lock_dir: './cache_data/locks'
    # where to store lockfiles for tile creation
    tile_lock_dir: './cache_data/tile_locks'
  image:
    resampling_method: bilinear # nearest/bilinear
    format: image/png

HTML - leaflet file code

 var crs = new L.Proj.CRS(
'OSGEO:41001','+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs',
{
resolutions: [8192, 4096, 2048, 1024, 512, 256, 128],
origin: [0, 0]})
var map = new L.Map('map',{
crs: crs
}).setView([21.46,78.3084],5);

L.tileLayer('http://localhost/mapproxy/tms/1.0.0/tms/osm_grid/{z}/{x}/{y}.png', {
maxZoom: 20,
minZoom: 0,
continuousWorld: true,
attribution: 'Map data © OpenStreetMap contributors, Imagery © Kartena'
}).addTo(map);
}
initmap();

Here is how the result looks like, against how it should be seen.

Result is looking like this

This is how result is expected and can be seen on MapProxy Demo


texts: Bivariate legend plugin throws NameError exception

texts:
While using QGIS 3.10 and 3.16 on Mac OS X Big Sur, the Bivariate legend plugin throws a NameError exception. How do I fix this or work around it?
The following stack trace is from QGIS 3.10. The 3.16 stack trace seems similar.
An error has occurred while executing Python code:
NameError: name 'Qgis' is not defined
Traceback (most recent call last):
  File "/Users/judithheher/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/BivariateLegend/bivariate_legend.py", line 400, in generate_image
    level=Qgis.Info
NameError: name 'Qgis' is not defined

Versions:
Python version: 3.7.3 (default, Mar 27 2019, 09:23:15) [Clang 10.0.1 (clang-1001.0.46.3)]
QGIS version: 3.10.11-A Coruña A Coruña, d2171173e4
Python Path:
/Applications/QGIS3.10.app/Contents/MacOS/../Resources/python
/Users/judithheher/Library/Application Support/QGIS/QGIS3/profiles/default/python
/Users/judithheher/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins
/Applications/QGIS3.10.app/Contents/MacOS/../Resources/python/plugins
/Applications/QGIS3.10.app/Contents/Resources/python
/Applications/QGIS3.10.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/site-packages
/Applications/QGIS3.10.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7
/Applications/QGIS3.10.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python37.zip
/Applications/QGIS3.10.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/lib-dynload
/Applications/QGIS3.10.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/site-packages/geos
/Users/judithheher/Library/Application Support/QGIS/QGIS3/profiles/default/python
/Users/judithheher/Documents/Home/SJSU/GEOG171/Week4


texts: How would you illustrate this in a report to a client?

texts:
The objectives are simple: They're adding a train station and the goal is to display how far you could get within 30 minutes before and after the new station is built. In this map that I'm working on here the circle size is calculated like this(all attributes with cost above 30 is removed, if there's time life when you get to a station we assume walking at 6 km/h for the remaining time and a limit at 20 minutes to prevent the circles from getting to big at origin.

Without the zoom layer the last circle is impossible to see since the size is really small. Any other good suggestions on how to illustrate this?

case when "cost" < 10
then
20*100*2
else
(30-"cost")*100*2
END

Green = Nothing's changed
Yellow = Discontinued public transport zone.
Blue = New public transport zone.

How would you solve this with only one map? Other solutions then mine wanted!


texts: Storing a geographic minimum spanning tree with PyShp

texts:
I have a list of geographic coordinates (longitudes, latitudes in decimal degrees) as follows:

-72.9833,10.6167
-63.1167,10.55
-63.1167,7.91667
-65.9667,10.3167
-67.4333,8.93333
-68.4,10.6167
-70.6667,9.53333
-72.8667,9.83333
-72.9833,10.6167
-60.5,3.93333
-56.4833,3.1
-78.9333,8.38333
-79.4333,9.15
-81.0167,7.66667
-81.1833,7.51667

and a list of edges connecting these points with a minimum spanning tree (computed using Prim's algorithm) as shown below:

0 0
5 1
4 5
11 5
10 11
14 10
3 14
13 3
15 13
8 4
7 8
12 13
2 12
9 2
6 9

I am trying to save the minimum spanning tree as a shapefile using the PyShp library (so that the MST can be displayed on a GIS).

Here is my code so far:

    import numpy as np
    from mpl_toolkits.basemap import Basemap
    import matplotlib.pyplot as plt
    import shapefile

    # Load coordinates
    coords = np.loadtxt('coords.dat')
    lon = coords[:,0]
    lat = coords[:,1]

    # Load graph edges
    edges = np.loadtxt('edges.dat')
    edges = np.array(edges, dtype=np.int32)

    # Create map
    map = Basemap()

    map.drawcoastlines()
    map.drawcountries()
    map.drawmapboundary()
    x,y = map(lon,lat)
    map.plot(x,y,'ro')
    plt.title('Minimum Spanning Tree')

    # Ths MST plotting algorithm works!
    edge_list = []
    n = len(edges)
    for i in range(n):
        t = edges[i,0]-1
        u = edges[i,1]-1
        x = [coords[t,0], coords[u,0]]
        y = [coords[t,1], coords[u,1]]
        #edge_list.append(x)
        #edge_list.append(y)
        #edge_list.append((x,y))
        edge_list.append([(x[0],y[0]), (x[1], y[1])])
        plt.plot(x,y,'-r')
    plt.show()

# Sort the list of edges
edge_list = list(sorted(edge_list))

# Save the minimum spanning tree
w = shapefile.Writer()
w.line(parts=[edge_list])
w.field("COMMON_ID", 'C')
w.record("Point")
w.save("mst")

The code above correctly plots the points and draws the minimum spanning tree connecting them, as shown in the figure below:

But the tree is not being correctly saved in the shapefile, as shown in the figure below (displayed by QGIS):

An alternative would be to have the tree saved as a KML file (compatible Google Earth), which could then be converted into a shapefile using the ogr2ogr tool of GDAL. I have not tried this as yet.

Any hints?

UPDATE 1: Following the suggestion of Marc Pfister (see comments below), I stored the coordinates as lists of tuples in edge_list: edge_list.append((x, y)). The results are shown below:

[([-61.049999999999997, -63.116666666699999], [10.4, 10.550000000000001]), ([-63.116666666699999, -63.116666666699999], [7.9166666666700003, 10.550000000000001]), ([-63.116666666699999, -65.966666666699993], [10.550000000000001, 10.3166666667]), ([-67.433333333299998, -65.966666666699993], [8.9333333333299993, 10.3166666667]), ([-67.433333333299998, -68.400000000000006], [8.9333333333299993, 10.6166666667]), ([-70.666666666699996, -68.400000000000006], [9.5333333333300008, 10.6166666667]), ([-70.666666666699996, -72.866666666699999], [9.5333333333300008, 9.8333333333299997]), ([-72.866666666699999, -72.983333333299996], [9.8333333333299997, 10.6166666667]), ([-63.116666666699999, -60.5], [7.9166666666700003, 3.9333333333299998]), ([-56.483333333300003, -60.5], [3.1000000000000001, 3.9333333333299998]), ([-78.933333333299998, -72.866666666699999], [8.3833333333300004, 9.8333333333299997]), ([-79.433333333299998, -78.933333333299998], [9.1500000000000004, 8.3833333333300004]), ([-79.433333333299998, -81.016666666700004], [9.1500000000000004, 7.6666666666700003]), ([-81.183333333299998, -81.016666666700004], [7.5166666666699999, 7.6666666666700003])]

The shapefile created by storing these data is still uglier than the previous one.

UPDATE 2: Fixed the storing of the (x,y) tuples in edge_list. It now works fine! Have updated the code above to include this fix.


texts: How to get the Lat/Lon of a  from a PBF file?

texts:
I'm working with a PBF file and I'm interested in getting each <way> tag alongside the coordinates of each subnode of the way.
Here's what information I'm working with and how I'd like to transform it:
<node id="31154978" version="2" timestamp="2011-09-16T08:36:33Z" lat="41.3444981" lon="9.2521174"/>
<node id="31155011" version="2" timestamp="2011-09-16T08:36:33Z" lat="41.3446635" lon="9.2524808"/>
<way id="4843086" version="9" timestamp="2020-02-18T11:00:33Z">
  <nd ref="31154978"/>
  <nd ref="31155011"/>
</way>

...

<way id="4843086" version="9" timestamp="2020-02-18T11:00:33Z">
  <node lat="41.3444981" lon="9.2521174"/>
  <node lat="41.3446635" lon="9.2524808"/>
</way>

I tried to do it programmatically where I read the file, then for each subnode of my way I try to get the correct <node> that matches the id and from there I can get the lat/lon.
It works more or less but it's super slow and uses a lot of RAM as I have to store each node to be able to later get back to them and fetch the coordinates.
I tried some external tools like osmium or osmosis to transform my PBF file to a geojson-compatible file but it seems to transform the <way> with the coordinates as the center of the road.
Is that something possible with pbf file? How can I get the full coordinates of a <way> easily?


texts: How to pass OAuth token to fetch web map without passing userid?

texts:
I'm using the ArcGIS JavaScript API 3.7 for web map service. The web map is based on Oauth2 authentication. I use java code to fetch the token. Then I need to fetch the map in javascript by calling esri.arcgis.utils.createMap(). But How to pass the token with this request?

I know ArcGIS JavaScript 3.7 offers an API called esri.IdentityManager.registerToken() that can pass the token will the request. But this method requires the user name as parameter. I don't want to put the username in the javascript code, as it's not safe.

Are there any other ways to pass the token with the request?

The api for esri.IdentityManager.registerToken() is here.


texts: Why can't I change precision of decimal number in new Field Calculator?

texts:
Can't click those fields. It worked in the past and now just don't. I can't figure it out.

And even it's set to three by default (and I can't change that), it ends up being infinite.


texts: Link non-spatial SQL Server query to a Geodatabase?

texts: Is it possible to build an object (of whatever kind) in ArcEngine code that links the geodatabase data with a non-spatial SQL Server query? How would you go about doing this?


texts:
I have ArratList<Segment> (N segments - array size is N) when:

public class Segment {
    Node vertex_1;
    Node vertex_2;
}

public class Node {
    double latitude;
    double longitude;
}

And one static point - origin point. How to work with GeoTools in order to get ArrayList<boolean> (size N) where each value is true/false regarding the question:

Is the 2 lines from origin point to both edges of the segment, intersects some other segments on their way?

For example here the result is {true, false, false} because the red line from origin point to second edge of segment 1 intersects on her way segment 3.

This question is similar to Algorithm to find intersections in a polar plane on StackOverflow but the difference is that here I want to work with GeoTools instead of implementing the algorithm which involving translation of Geographic measurement units (lat/lon) into polar plane and performing some math calculation such as cross products - Not difficult but has potential for bugs, and if there is already ready open source library it's better to use it.

Because this question involving code solution, it asked also as Calculating intersection between lines with GeoTools (java) at StackOverflow.


texts: Python Syntax to ListTables if does not contain *

texts:
I am trying to create a list of tables on the condition that they do not contain 1 of 2 strings which would indicate the table is a view.  I have successfully built the table list with only 1 exclusion, but I can't figure out the syntax to include both exclusions.

The following line correctly builds a list excluding tables that contain the '.VW_' string.

tableList = [tb for tb in arcpy.ListTables('*.' + u + '.*') if '.VW_' not in tb]

And the following line correctly builds a list excluding tables that contain the '.TEMP_' string.

tableList = [tb for tb in arcpy.ListTables('*.' + u + '.*') if '.TEMP_' not in tb]

I have tried to combine the logic to exclude a table if it contains '.VW_' or '.TEMP_', but I can not seem to the syntax right.  The result of the following line includes all tables even if they have '.VW_' or '.TEMP_' in the name.

tableList = [tb for tb in arcpy.ListTables('*.' + u + '.*') if '.VW_' not in tb or '.TEMP_' not in tb]


texts: How to change the Help Description in Modeler?

texts:
I´ve created some Models in QGIS 2.18 and the help in the Model shows only what is in the description field. Thats the way I want it.

But I migrated to QGIS 3.4, and now when I type something in the Modeler help, all descriptions like input, output and so on, are displayed.

Is there a way to limit the information displayed, to be like in 2.18?


texts: Calculate areas of overlapping polygons from two shapefiles

texts:
In Why does intersect not calculate polygon areas? the answer is "The area values in attribute table are just random attributes in the eyes of QGIS. It is therefore perfectly normal that these attribute values are not manipulated when running Intersect.

Intersect is provided by ftools plugin (a Python plugin). You can expand the tool's code to calculate the area of the resulting intersecting areas."

What would be the "code to calculate the intersections' areas" and to export them? It would be ingenious to do all tasks in ONE tool.


texts: Distance Matrix with 2 point shapefiles and one street network

texts:
I have two point shapefiles (start/destination) and one sort of street network shapefile.
I would like a to create a distance matrix from the start-shapefile to the destination points, using the street network file.
It does not necessarily have to be "shortest path", because mostly there is only one way anyway.

I am using QGIS 3.21 at the moment, but I have 2.18 also, so I don't really care on which one.

Is there a proper solution?


texts:
I've started learning pyQGIS and the first step was to work on easy things like to add raster and vector in QGIS. So, I found this code in Loading raster layer using PyQGIS? and it works fine and is easy.

from qgis.core import QgsRasterLayer
from PyQt4.QtCore import QFileInfo

def StringToRaster(raster):
    # Check if string is provided

    fileInfo = QFileInfo(raster)
    path = fileInfo.filePath()
    baseName = fileInfo.baseName()

    layer = QgsRasterLayer(path, baseName)
    QgsMapLayerRegistry.instance().addMapLayer(layer)

    if layer.isValid() is True:
        print "Layer was loaded successfully!"

    else:
        print "Unable to read basename and file path - Your string is probably invalid"

raster = 'C:/home/user/Desktop/output.tif'
StringToRaster(raster)

However, I have some question : this code works with one absolute path to one raster but, if it exists a list of paths like this:

raster = ['C:/home/user/Desktop/output.tif','C:/home/user/Desktop/output2.tif','C:/home/user/Desktop/output3.tif','C:/home/user/Desktop/output4.tif']

How to make this code works such it adds all the rasters in the list?


texts: How do I move and/or rename a file geodatabase?

texts:
I am expanding the scope of an existing file geodatabase so I would like to move it up a folder hierarchy as well as rename it.

For both of these operations, is it as simple as moving the folder via Windows explorer as well as renaming the folder via Windows explorer?

Complications I'm trying to avoid are if there is some metadata that will be lost within the gdb if the name changes. I can restablish links in mxd files that reference the gdb easily of course.


texts: Scoring fields between 0 and 1 using ArcGIS field calculator?

texts:
I am starting to get to know ArcGIS 10.2.1 and I need to the following, but I am not sure how to build the expression into the FIELD CALCULATOR:

I need to add first a new field called Near_distance_to_LAriver which it should score between 0 and 1  the nearst distance of 4 sites to Los Angeles River, using the default field NEARDIST values, using the field calculator, I need to score them between 0 and 1 Being 1 the nearest site to the river (this means: the shortest distance). This new field should be the tipe FLOAT, but the to score it I dont know which one should be the expression at the Field calculator


texts: Understanding measurement units for MODIS data?

texts:
I am working on MODIS data in ESA SNAP, specifically the Terra sensor, and I'm having a hard time understanding the unit of measurement of the data provided.

I have downloaded Level 1B data (Calibrated Radiances, 1km resolution). After loading the product, I can easily open the image from any band. Once opened, if I open up Menu/View/Tool Windows/Pixel Info, besides geographic information, what I see is the "Band" value.

What is this value?

It is not expressed in any unit of measurement. Now, the MODIS website suggests that they are radiances (and this is what I would expect, as the name of the product suggests). On the other hand, if I open the "Product user guide" in the same page, I would expect to see these famous "Scaled integers", as briefly explained in page 31, section 5.4 of the guide.

The following is an example of what I mean. In the "bands" field, I would expect integers (the fields are two because I tiled two images). Also, tie-point grids are expressed in no unit. Is there an SI unit that I'm not aware about?

I have appealed to people who know a lot more about me about SNAP and MODIS data, and, with @chryss' help, I've come to the answer. I looks as if SNAP automatically processes the "Scaled integers" to output reflectance (dimensionless, %) for the reflective bands and calibrated radiances (in W/(m^2*ster*μm)) for the infrared.


texts: Labels missing after conversion from *.mxd to *.kml?

texts:
I've been trying to convert a map from .mxd to .kml for google earth. The map has municipalities labeled as well as train stations. When I save it and open it in Google Earth, only the train station labels appear.

I have tried moving the municipal layer up in the table of contents before saving. I have also tried taking out the rail stations labels and only having the municipal labels. Neither has done the trick.


texts: Infowindow close on featureOut /hoverOut

texts:
I have been working on getting infowindows on my cartodb layers, but i am not sure what I am doing wrong as I want the infowindows to automatically close on feature out and only pop up when featurein or on hover regardless of the sublayer. A simple scaled down approach is as shown here. http://bl.ocks.org/jmwenda/b472b0a932fc5b895a5e59dadac18192. I have tried other methods such as shown with very little success. Any pointers on what I going on here or how to trigger the infowindow to automatically close?

  cartodb.createLayer(map, markersLayerSource)
    .addTo(map)
    .done(function(layer) {
             var publishedLayer = layer.getSubLayer(0);
             publishedLayer.setInteraction(true);
             var infowindow = publishedLayer.infowindow;
             cdb.vis.Vis.addInfowindow(map, publishedLayer,['cartodb_id','cdb_list'],{triggerEvent: 'featureOver'});

             publishedLayer.on('featureOver', function (e, latlng, pos, data, layer) {
                      //$(document).bind('mousemove', event);
                      $('.cartodb-popup').show();
                      //infowindow.set('visibility', true);

             });

             publishedLayer.on('featureOut', function (e, latlng, pos, data, layer) {
                      $('.cartodb-popup').hide();
                      //$(document).unbind('mousemove', event, false);
                      //infowindow.set('visibility', false);
             });

    });


texts: Adding name of used coordinate system using alabel in QGIS Layout?

texts:
Is there a way to display the name of the used CRS in the layout?

The expression "project_crs" displays the code of the used CRS but I'd like to use the common name of the reference systems.

This question was asked and answered before as Adding coordinate system to layout in QGIS composer?, however the api the answer used changed with QGIS 3.0 and the solution is not valid anymore.


texts:
I used to be able to Georeference in QGIS but now I can't. These are the steps I go through. At the end of it nothing happens. No layer is added to my Layers Panel and no document is created in the nominated directory. Can anyone see what I am missing or doing wrong?

Click on ‘Web’ (in the top menu), then ‘OpenLayersPlugin’
Zoom into NE coast of Australia and find 5 points on the OpenStreetMap that I can match up on the scanned map.
Click on ‘Raster’ in the top menu and then ‘Georeferencer’
Click on the ‘Open Raster’ icon
Make sure that the CRS of the scanned map is EPSG:3857 (ie. same as OpenStreetMap)
Click a point on the scanned map, then choose ‘From map canvas’, then click the corresponding point on the base map (OpenStreetMap), and click OK
Repeat the process described at #6 five times (ie. create 5 registration points)
Click on ‘Transformation Settings’ and use following settings:

         Transformation type: ‘Thin Plate Spline’
         Resampling method: ‘Nearest Neighbour’
         Target SRS: EPSG:3857 WGS 84 Pseudo Mercator
         ‘Load in QGIS when done’ = ticked


texts: ogr2ogr in Windows cmd prompt unable to find driver 'PostgreSQL'

texts:
I can successfully execute the following command in the OSGeo4W shell:

ogr2ogr -f "PostgreSQL" -a_srs "EPSG:102733" PG:"host=localhost user='postgres' dbname='nameofdb' password='password' port=5433" "Water and Sewer Map.dxf" -nln eastover

but whenever I open up a blank Windows CMD and navigate to the folder where my file is and run the same command, I get this:

Unable to find driver `PostgreSQL'.

But when I run this in the Windows CMD:

ogr2ogr --version

I get this:

GDAL 2.0.1, released 2015/09/15

My goal is to have this as a batch file so I can run it from the native CMD prompt. Why is it saying the I don't have the PostgreSQL driver on the Windows CMD and how to I fix it?

UPDATE

Per AlecZ's recommendation, I have taken the o4w_env.bat script and appended my command to the end to look like this:

REM Make parent of this script location our current directory,
REM converting UNC path to drive letter if needed

pushd %~dp0
cd ..

REM set OSGEO4W_ROOT to short path version
set OSGEO4W_ROOT="C:\OSGeo4W64"

REM start with clean path
set path=%OSGEO4W_ROOT%\bin;%WINDIR%\system32;%WINDIR%;%WINDIR%\WBem;%OSGEO4W_ROOT%\share

for %%f in ("%OSGEO4W_ROOT%\etc\ini\*.bat") do call "%%f"

popd

cd network share:\folder\folder\folder

ogr2ogr -f "PostgreSQL" -a_srs "EPSG:102733" PG:"host=localhost user='username' dbname='databasename' password='password' port=5433" "Water and Sewer Map.dxf" -nln eastover

pause

But now I am getting this error:

ERROR 4: Unable to open EPSG support file gcs.csv.
Try setting the GDAL_DATA environment variable to point to the
directory containing EPSG csv files.
ERROR 1: Failed to process SRS definition: EPSG:102733

I have updated the environment variable as suggested in the error and added ;%OSGEO4W_ROOT%\share to the end of the clean path setting, but neither seems to be working. Any other suggestions?


texts: Duplicating Features To Create A Pattern

texts:
I'm currently trying to take the 'Shifting Features' script from ArcPy Cafe and this SE response and apply it to a similar script that will one by one, through all the features in a feature class, copy the first feature to another feature class, then create a copy of that feature X times with a set interval and angle.

So with this I should be able to take, say a feature class containing two lines, and create 5 lines, 3 feet apart at a 45 degree angle in a new feature class.

I figured this should be easy by nesting cursors, one for searching the original FC, one for inserting rows in the new FC, but I apparently don't know much about cursors in general, or SearchCursor or InsertCursor in general. The InsertCursor command runs the correct number of times, and creates the correct number of features, but the features have no geometry.

Can anyone help me figure out what I might be doing wrong? Thanks!

import arcpy, math

# all of these will eventually be made arguments

fc = r"c:/geometry/polylines.shp" # source feature class
fc2 = r"c:/geometry/polylines_shift.shp" # destination feature class
c_shift=3 # desired overall offset (hypotenuse)
no_features=5 # number of features to duplicate between offsets
angle=45 # angle of offsets

# given the angle and the hypotenuse, determine the x and y coordinate shift
x_shift=c_shift * math.cos(angle)
y_shift=c_shift * math.sin(angle)

# I know this will copy the existing features from fc to fc2
#This will eventually be replaced with a command to create a blank feature class once I get the meat of the script fixed
if arcpy.Exists(fc2):
    arcpy.Delete_management(fc2)
arcpy.Copy_management(fc,fc2)

i_cursor = arcpy.da.InsertCursor(fc2,("SHAPE@XY"))
with arcpy.da.SearchCursor(fc, ["SHAPE@XY"]) as cursor:
    for row in cursor:
        for n_feat in xrange(1, no_features+1):
            i_cursor.insertRow([[row[0][0] + (x_shift*n_feat or 0),
                                 row[0][1] + (y_shift*n_feat or 0)]])
del row, cursor, i_cursor


texts: Prevent GeoPackage from saving when it is not edited in QGIS

texts: So in QGIS, GeoPackage will be automatically saved when it is read even if you didn't edit. How to avoid this from happening?


texts: Building zig-zag lines (Drone flight path)

texts:
I need to create a line path like the one below.
I have to create a flight plan for an inspection with drone of a solar plant and I know the length of every single row and the distance between row, so I need a fast and automatic method to create the flight path.

There's any way to do that in QGIS using plugin or other function?


texts: Using Python script tool in ModelBuilder?

texts:
i have created the following tool to add two fields in a table and do a calculation in order to put the tool to ModelBuilder.

For some reason when i am going to put it as an input in an another tool it doesn't appear.

I think i have done something wrong with
 1. SetParameter and also
 2. i haven't managed to save it as a feature class and it's only in memory.In addition i don't have a gdb or an mxd.

Also i set set parameter as an output feature class at tool's parameters but it didn't make any difference.

I am using ArcGIS 10.2.2, any suggestion?

 # Import system modules
import arcpy
from arcpy import env
arcpy.overwriteoutput = True

# Set environment settings
inWorkspace = arcpy.GetParameterAsText(0)
inFeature = arcpy.GetParameterAsText(1)
outFc = arcpy.SetParameterAsText(2, "")
env.workspace = inWorkspace

# Set local variables
fieldName1 = "bear_60"
fieldName13 = "distance"

# Execute AddField for 2 new fields
arcpy.AddField_management(inFeature, fieldName1, "DOUBLE", "", "", "", "", "NON_NULLABLE")
arcpy.AddField_management(inFeature, fieldName13, "DOUBLE", "", "", "", "", "NON_NULLABLE")

# Calculate Field bear_60
inTable = inFeature
fieldname = fieldName1
expression = "getCalc(!aspect!)"
codeblock = """def getCalc(aspect):
    if aspect > 300:
        bear_60 = (aspect + 60) - 360
    else:
        bear_60 = aspect + 60
    return bear_60"""
# Execute CalculateField
arcpy.CalculateField_management(inTable, fieldname, expression, "PYTHON_9.3", codeblock)

# Calculate Field distance
inTable = inFeature
fieldname = fieldName13
expression = "getCalc(!distance!)"
codeblock = """def getCalc(distance):
    return 10000"""
# Execute CalculateField
arcpy.CalculateField_management(inTable, fieldname, expression, "PYTHON_9.3", codeblock


texts:
I can use ArcGIS Desktop model builder to leverage the ArcPad Check In Tool.  When I run the model it executes successfully.  When I export the model as python I get this code:

import arcpy

# Local variables:
Test02_axf = "C:\\aWorkspace\\ESA\\ArcPad\\AXFs\\Test02.axf"

# Process: ArcPad Check In
arcpy.ArcPadCheckin_ArcPad(Test02_axf, "", "")

However, when I run this code in the python window I get:

Runtime error
Traceback (most recent call last):
File "<string>", line 8, in <module>
AttributeError: 'module' object has no attribute 'ArcPadCheckin_ArcPad'

And yes indeed arcpy has no function named ArcPadCheckin_ArcPad.  So model builder must be exporting erroneous code.

Does someone out there know what the issue is here?  Is there really anyway to do this operation with python?


texts: How to redesign confusion matrices from Ecognition with python and seaborn?

texts:
I am using eCognition manual classification as a baseline for performance comparison with other 3D classification methods. For each method, I am creating confusion matrices in python with sklearn and seaborn. I am trying to transform the confusion matrices resulting from the eCognition accuracy assessment (see image 1) into seaborn generated confusion matrices so they all have the same pattern for publication i.e with the same seaborn colour palette as in image 2.

The code I used to create the confusion matrices for the other methods is below:
Here, the y_true is the annotated data (ground-truthed) and the y_pred is the result from the classification. The input is a csv file where each pixel has its manual annotation (y_true) value and its y_pred value.
def plot_cmatrix(y_true, y_pred, figsize=(10,10)):
    data = confusion_matrix(y_true, y_pred)
    df_cm = pd.DataFrame(data, columns=np.unique(y_true), index = np.unique(y_true))
    df_cm.index.name = 'Actual class'
    df_cm.columns.name = 'Predicted class'
    plt.figure(figsize = (10,7))
    # plt.title("Model F \nAccuracy: {0:.3f}".format(accuracy_score(y_true,y_pred)))
    plt.title("Model F")
    sns.set(font_scale=1.4)
    sns.heatmap(df_cm, cmap="Blues", annot=True, fmt='d', annot_kws={"size":16})
    tn, fp, fn, tp = confusion_matrix(y_true, y_pred).ravel()
    print("True negative: ", tn,"False Positive: ", fp, "False Negative:", fn,"True Positive: ", tp)

So I want to apply the function with eCognition classification results, but I am not sure how to get information on y_pred and y_true information there. This brings me to the three questions below:
Would anyone know how can I get these values from the classification in eCognition?
Any suggestions are welcomed
If you need any further information please feel free to contact me


texts: Maplex not placing labels in area that appears to be available?

texts:
I am trying to get ArcMap 10.2.2 to place my labels where I want them using the Maplex label engine. I want to stop the labels from both overlapping other features in the same layer, and overlapping a text annotation layer. The result I am getting with my current settings is this:

As you can see, the two red labels overlap the annotation text (and a bit of a feature too) even though at least one of them should fit into the white available area. It seems like Maplex doesn't consider the low priority user defined zones.

I have the following settings:

Position: Offset, and user defined zones where to place them. (Starting from 12 o'clock the wheel reads 7, 3, 1, 4, 8, 5, 2, 6.)
Label offset of 1 mm, maximum 500% of desired.
No fitting strategy in use.
Never remove as conflict resolution. If I uncheck this, the two red labels are not shown at all.
0% label buffer.
The feature class has a weight of 1000 (I have tried various other values, but it does not help).
The text annotations has weight "high".
Using "Best", not "Fast".

Please let me know if any other settings are of interest, and I will add them.

I know that label fitting is more of an art than a science, but I want MapLex to utilize all the available space.


texts: Creation of a raster dataset from polylines

texts:
I face an issue with the Raster tool. I have Polylines that I want to convert to Raster. The spatial reference coordinate system is World Mercator (WGS 84).

I used the "Convert Polyline to Raster" tool. But the problem is that I expected to get an attribute table with all the cells (I defined the cell size as 1000) with values corresponding to the ID of the line crossing the cell... However, I get an attribute table with all values greater than zero and therefore not all the cells.

Is this was linked to the Coordinate System, and if so, do I need to re-define the project and use a system like Lambert Conic (Projection System)?


texts: How to display 33k placemarks with custom symbols/labels?

texts:
I want to overlay address numbers on our embedded Google Map (V3 api).   Even though the KML is under 2 MB, it exceeds the 1,000 feature limit; there are over 33,000 placemarks.

Fusion tables solve the size restrictions, but do not support custom symbols.  Clicking on a marker to display the address number in a bubble is not acceptable.  Did I miss something obvious that supports creating labels?  I can see the address numbers in the Name column...

Are custom image overlays the answer?
http://oregonstate.edu/campusmap/
http://www.mail-archive.com/google-maps-js-api-v3@googlegroups.com/msg09037.html

Please share tips on how to create the png image from shapefile vector points.  I have access to ArcGIS 10 if that's an option.

Thanks!


texts: I maintain metadata as I always have....inconsistently and sporadically.. I know...booh to me. Anyways, I know python can pull all kinds of info out of mxd in  V10. Has anybody come up with some scripts to populate basic stuff like layers and paths to mxd metadata? Seems like this should be pretty easy, but I would really like to know if there is any reason this won't work and if anybody has written anything along these lines to get me started in the right direction.


texts: Creating single Geoserver layer from two contiguous JPG+WorldImage files?

texts:
I have two sets of files for two contiguous areas. I've been asked to merge these two areas into one and serve the resulting data set in Geoserver. What I have in each directory is what I reckon is a WorldImage data set generated with ArcGIS: I have a JPEG file, and others with the suffixes .jgw, .prj, .rrd, .aux, and .aux.xml.

I was able to convert each to a GeoTIFF, then merge the resulting GeoTIFFs with gdal_merge.py, with more quality loss than the boss likes, and I was told it should be possible to avoid going through merging and just tell Geoserver to create a layer covering both involved and contiguous areas.

I'm not sure that's possible, but I could be wrong. Is that the case?


texts: Openlayers - Geoserver - tranforming srs of layers to match the base layer

texts:
I have the following two layers and Open Street Map (OSM) layer. A part of the script is given below. My problem is when I view the results from the browser both the layers fall on a wrong coordinates. I think this is because they are originally ESPG 3006 but the base layer (OSM) is ESPG:3857. Now I want to transform the two layers to match with the base layer ESPG. Can i do this without using any GIS software? Please Help.

vagar_alla = new OpenLayers.Layer.WMS(
                "vagar_alla", "http://localhost:8000/geoserver/myProject/wms",
                {
                    height: '489',
                    width: '512',
                    layers: 'myProject:vagar_alla',
                    styles: '',
                    srs: 'ESPG:3006',
                    format: 'image/png',
                    transparent:true,
                    tiled: 'true',
                    tilesOrigin : map.maxExtent.left + ',' + map.maxExtent.bottom
                },
                {
                    buffer: 0,
                    displayOutsideMaxExtent: true

                }
            );
adresspt = new OpenLayers.Layer.WMS(
    "adresspt", "http://localhost:8000/geoserver/myProject/wms",
    {
        height: '512',
        width: '458',
        layers: 'myProject:adresspt',
        styles: '',
        srs: 'ESPG:3006',
        format: 'image/png',
        tiled: 'true',
        transparent:true,
        tilesOrigin : map.maxExtent.left + ',' + map.maxExtent.bottom
    },
    {
        buffer: 0,
        displayOutsideMaxExtent: true
    }
var osm = new OpenLayers.Layer.OSM();

    map.addLayers([osm,vagar_alla,adresspt]);
);


texts: How to get unique values and cell counts of a raster in QGIS 2.6?

texts:
I have a soil layer that I have clipped to a polygon area and I want to obtain the unique cell values (eg. each soil type) and the cell count for each value.

I have a lot of experience with ArcMap (where this is a simple procedure) however I want to learn to use QGIS. I have looked online and can't find an answer. Is there an easy way of doing this?


texts: Would any of these show the depth of lakes in great detail? i.e. a LiDAR at 25cm, will it show the detail of the lake bed?


texts: Load ArcGIS Server layers into web using ArcGISDynamicMapServiceLayer JavaScript API

texts:
Don't know what I'm doing wrong; trying to load layers from ArcGIS Server  into web using 'ArcGISDynamicMapServiceLayer', can see the basemap but not the layer.

see code below:

<html>

<head>

<
<title>FeatureLayer</title>

<link rel="stylesheet" href="https://js.arcgis.com/3.19/esri/css/esri.css">

<script src="https://js.arcgis.com/3.19/"></script>

<script>
require([
    "esri/map",
    "esri/layers/FeatureLayer",
    "dojo/domReady!"
  ],

function(
    Map,
    FeatureLayer
  )
 {

 // Add data from ArcGIS Server

var featureLayer = new esri.layers.ArcGISDynamicMapServiceLayer("https://herse.chelmsford.chelmsfordbc:6443/arcgis/rest/services//Boundaries/MapServer/2");

map.add(featureLayer);

  });
</script>


texts: Is it possible to copy and paste features in an online Edit function?

texts: Using the ESRI JS API, is it possible to select a feature and then create buttons to copy and paste that feature?  This would be done using the Editor Widget.  If its possible to do this, but only retain the shape (loosing the attributes) I could live with that.  If its not possible using the ESRI JS API, can it be done simply with Dojo?


texts: How to group adjacent polygons with similar orientations?

texts: I have a layer of about 40k polygons, each representing an individual agricultural field.  I want to group adjacent fields (within 20meters) that are within 15degrees of orientation of each other.  I need to be able to measure how many of the 40k fields can be placed into groups with other fields that are contiguous and parallel.   I have already calculated the orientation of each field and converted them to a range of 0-180.  I am fairly familiar with ArcMap, however, I am just recently learning to apply my limited coding skills.


texts: Merging/Exporting Data Based on Location ArcGIS 10.3

texts:
I am trying to identify demographic characteristics of areas within a 1-mile radius of schools.

I have one shapefile with block level census data and another with the school locations.  The files do not have a shared field.

Is there a way to merge the data to create a shapefile that contains the census blocks (and accompanying data) and the schools?  There would obviously be repetition among the schools, since multiple census blocks fall within each school's radius. There would also be some overlap/repetition of census blocks, since some schools are close enough to each other to make some radii overlap.

I am pretty sure there is a simple solution that I learned at some point, but my GIS skills are a bit rusty and I'm having difficulty.


texts: Polygons in one layer have different attributes

texts:
Is it possible in PyQGIS to have in one layer more polygons and some would have one set of attributes and another different?

E.g. Polygon_1(name, size, average_age), Polygon_2(name, size, average_height), Polygon3(same as Polygon1)...

I imagine that I would create two tables and when I create a new polygon I somehow link it together.
Is anything like that possible? Or do I have to have one layer for each set of attributes?


texts: How to combine GEOMs or linestrings from different rows

texts:
I am working with a set of data that includes common IDs and overlapping linestrings/geoms. I would like to combine the overlapping lines based on their common IDs, preferably in postgres. I can run the combo process on either format. The linestrings are lat/long.

The data as linestring (lat/long) looks like this (with each # being a lat/long pair):

ID     geometry-points
11     1 3 56
11     12 45 7
11      234 6 8
12       3 5
12     1
12           78

And would return like this

11     12345678
12     13578


texts: ArcGIS Enterprise SSL error when loading hosted feature service

texts:
Long story short: An item from our portal server is skipping around the portal URL, and loading as its own resource. This is causing a problem because the server does not have its own certificate. It should be bundled up with the portal resources.

Note: https://internal.gis.bisnet.com is on https://server01.ina.bisnet.com

We have ArcGIS Enterprise deployed. Everything is secure, however when a hosted feature service is loaded in a new web map, the site is flagged as insecure.

In Google Chrome Developer tools, I loaded the site to understand what insecure resource was being loaded.

It notes that my portal site (https://internal.gis.bisnet.com --example) IS secure.
However, it is also showing the name of the server that portal is on (https://server01.ina.bisnet.com) as an INSECURE resource. This is interesting, because everything should be coming through the portal URL (internal.gis.bisnet.com).

I took a look at the Google Chrome sources list, to determine what resource was loading from server01 (skipping around portal).

It is an item from our Webadaptor (called team)
It appears to be something under the "self" properties

How do I modify my settings so this resource travels through portal, and does not appear as its own resource?


texts: How to read Landsat images from Amazon S3 using Python?

texts: I would like to read and process Landsat images from Amazon S3 using Python


texts: Python for ArcGIS: Code refuses to accept a defined variable as input

texts:
I'm learning to use arcpy, and I am starting to try and script some geoprocessing tools.

While trying to use arcpy.Select_analysis, I found that the code does not accept a defined variable as input. If (and only if) I directly type my file location into the initial input, will the code work. Once I try to define the file location as a variable, the script tells me that the dataset is not supported. This doesn't make sense to me - I copied the file location correctly, nothing has changed except that I defined it as a variable. I even changed the variable name a number of times, just in case.

This is the initial code:

import arcpy
inLayer='file_path'
outLayer='output_dir'
arcpy.Select_analysis('inLayer', 'outLayer', '"CNTRY_NAME"=\'Mexico\'')

The only way I've managed to make it work is to do:

import arcpy
arcpy.Select_analysis('file_path', 'output_dir', '"CNTRY_NAME"=\'Mexico\'')

Edit:
Actually I was using both inLayer and outLayer with arcpy.GetParameterAsText()
But to try and figure this out I changed it all to filenames instead.
My paths were file_path=r'C:\LabData'

Again, it worked when I put it directly into the code, but not when I defined it as a variable.

This error carries over when I define the variable to accept user input; no matter what I do, the initial step refuses to take any input and just says that the dataset is not supported.

I couldn't find anything similar to that (or I was searching the wrong terms) on here so I was hoping someone could help.


texts:
I'm trying to create a new vector layer and store it in a MySql-database, but when I test the connection I get an error with no detailed information at all. If I give an erroneus password I get a message that says:

  MySQL connect failed for:
  db_qgisgeodata,host=localhost,port=3306,user=xxxx,password=ppppppp
  Access denied for user 'uuuu'@'localhost' (using password: YES)

Because I get a detailed error message if i give the wrong login credentials I think QGIS is able to connect to the database, but there is something more to it.

The user has full permissions on the database. The problem is the same if I try with the root-user.

The user has full permissions on the database. The problem is the same if I try with the root-user.
The database is empty.
Windows 7x64
QGIS 1.8.0
MySql 5.5.28
it's no problem to connect with the database from other applications.

I guess I'm missing something quite obvious here, but I can't figure out what. Any ideas anyone?


texts:
I'm using PostGIS with OSM table of buildings.

How can I compute isovist from single location point while only the buildings (represented as geom's) are the obstacles?


texts: Calculate mean value for each raster and convert output array to data frame

texts:
How do I calculate mean value for each raster and convert output array to data frame?

r <- as.data.frame(cellStats(x,mean))


texts: How can I get a web-hosted kmz file to project correctly when added as a layer to an arcgis online map?

texts:
I have a kmz file hosted on dropbox that I would like to display in a map through arcgis online. The kmz file was generated from an esri grid raster. I have tried altering the raster projection prior to converting it to a kmz but this appears to have no effect.

You can see the arcgis online map here.

In google maps, the kmz is displayed without any errors:

Has anyone encountered this problem before and perhaps encountered a solution as well? I would just use google maps to display this data instead of arcgis online but it seems they won't be supporting kml files after February.


texts: Add New Project -> Class Library (Engine) targets the .NET 2.0 framework

texts:
Any one have comments on why the ArgGis Extensions in the 10 framework still target the .NET 2.0 framework?

For clarification do the following steps in Visual Studio:
Add New Project -> Class Library (Engine)  -> Click properties -> see that it targets the .NET 2.0 framework.


texts: How to ignore invalid geometries in performing spatial join in QGIS 3.0?

texts:
I am attempting to perform a very simple join attributes by location between two shapefiles in the same projection. When I run it, I get an error message indicating that the algorithm could not be run because of errors in geometry. It further states that I can ignore errors by choosing that option. But I don't see that option in the window to run spatial join.

Error message [running a French QGIS here]:

  L'entité (1771) a une géométrie non valide. Veuillez réparer la
  géométrie ou modifier le paramètre de traitement dans l'option
  "Ignorer les fonctionnalités d'entrée non valides". Execution failed
  after 2.58 seconds

Screenshot of options for running join by location:

One of the shapefiles (the target shapefile) has features where one polygon is located inside the other polygon. (When I run validate geometry, the error messages look like this: Le polygone 2 est à l'intérieur du polygone 0 [Polygon 2 is inside polygon 0].) This isn't an geographical "error" in the sense that a single feature in the database has some smaller polygons within larger polygons that should be considered a single feature and should not be considered as two separate features in the data table. Thus, I can't "fix" the geometry error without creating a database error. As mentioned in this post, this is common in my file and difficult to work around.


texts: Openlayers WFS Vector Layer with Popups: how to setCenter on point when you click it?

texts:
Working code I have so far is below. Lots of airports on the map, you click an airport icon, the popup appears. No problems there.

Goal: to click a popup, and based on the lat/long data of that particular point, the map will center on that clicked point.

wfs.events.on({
    featureselected: function(event) {
            var feature = event.feature;
            feature.popup = new OpenLayers.Popup.FramedCloud("box",
                    feature.geometry.getBounds().getCenterLonLat(),
                    null,
                    '
Airport: '+feature.attributes.Name +'
 ICAO: '+feature.attributes.ICAO+'
 IATA: '+feature.attributes.IATA+'
',
                    null,
                    true
            );
            while( map.popups.length ) {
                    map.removePopup( map.popups[0] );
            }
            map.addPopup(feature.popup);
            }
    });

Not sure how to go about this after browsing through other questions/looking at documentation. Also:

feature.geometry.getBounds().getCenterLonLat()

What exactly is the point of this? Based on the docs it's for {OpenLayers.LonLat} but it's not like it is reading the data from there?

Any input appreciated.


texts: Finding the middle point of a line in QGIS - Where to place input file name

texts:
I am trying to follow the python code below (Finding middle point (midpoint) of line in QGIS?) to create midpoints of line segments in QGIS 1.8. I am receiving errors executing the code. Unfortunately, I have no background in Python.

>>> def pair(list):
...    '''Iterate over pairs in a list '''
...     for i in range(1, len(list)):
  File "<input>", line 3
**for i in range(1, len(list)):**
^
IndentationError: unexpected indent**

another error:

>>>         yield list[i-1], list[i]
**File "<input>", line 1
yield list[i-1], list[i]
^
IndentationError: unexpected indent**

and finally

**Traceback (most recent call last):
File "<input>", line 1, in <module>
NameError: name 'mylayer' is not defined**

Pyhton code:

def mid(pt1, pt2):
   x = (pt1.x() + pt2.x())/2
   y = (pt1.y() + pt2.y())/2
   return QgsPoint(x,y)

def pair(list):
    '''Iterate over pairs in a list '''
    for i in range(1, len(list)):
        yield list[i-1], list[i]

def create_geometry(point,pr):
    # create geometry record
    seg = QgsFeature()
    seg.setGeometry(QgsGeometry.fromPoint(point))
    pr.addFeatures( [seg] )

# memory layer
pt_layer = QgsVectorLayer("Point", "midpoint", "memory")
pr = pt_layer.dataProvider()

for elem in mylayer.selectedFeatures():
    line = elem.geometry()
    for seg_start, seg_end in pair(line.asPolyline()):
        line_start = QgsPoint(seg_start)
        line_end = QgsPoint(seg_end)
        # midpoint
        midpt= mid(line_start, line_end)
        # add midpoint point to layer
        create_geometry(midpt,pr)
        pt_layer.updateExtents()

QgsMapLayerRegistry.instance().addMapLayers([pt_layer])


texts: ArcGIS WebApp Near Me widget is missing features

texts:
When I set up the Near Me widget in my Web App configured to 'Only draw selected layer results' the map view does not display all of the features in my layers (even before I have selected a layer in the results list). I would like to use this option but I cannot trust the results.
Have a missed something?
Everything displays fine when this option is not selected or when not using Near Me.


texts:
I have some big mosaic files (about 3GB, 50k x 50k) and use gdal_retile.py to generate the pyramid tiles for geoserver. It takes very long time to generate all tiles (about one day in my computer).

After googling, I didn't find a parallel version of gdal_retile.py. Are there other methods to generate tiles for big mosaic with higher speed?


texts: How to remove unwanted parts of a bigger map that is outside the smaller map above it?

texts:
This is what I'm working on:

How can I remove the parts of the land cover map that is outside the Baao map?

I can't use the filter function since the land cover map does not have location data:


texts: How do I solve this dependency error in openSUSE 13.1 when installing QGIS (master)?

texts:
I am running OpenSUSE 13.1 and downloading via Yast I get the following error:
nothing provides libpq.so.5()(64bit) needed by qgis-master140102-1.2.x86_64

I get three of these errors, for each of the three qgis downloads (master, grass plugin, master-devel).

I have two options under Conflict Resolution:
1-do not install...
2-break...by ignoring some of its dependencies

What would you folks recommend? Is there something wrong with QGIS in openSUSE?


texts: Deleting Custom Label Placements in ArcGIS Pro

texts:
I was doing some trial and error by creating Custom Label styles and placements and have achieved what I set out for. Problem is, I've got a few 'bad' custom label placements I want to delete and only keep the good ones.
I know how to delete Styles where do I delete the custom label placement icons?


texts:
I have two layers, first one is the ROAD layer and second one is the PARCEL layer (both polygon layers) as shown in Figure 1. I can get the data both in DXF and shapefile formats.

My task is to compute the area of intersection between ROAD and PARCEL layer, this is the easy part. I can compute intersection using QGIS or Geopandas easily. However, the difficult part is creating maps for each of the parcel. Sometimes, I have to create more then a hundred maps for each project. For mapping, there is a template that I have to use, which is similiar to the one shown in Figure 2. Also, some attribute data should be included in the map, such as owner of the parcel.

These maps should be in both PDF and DXF formats. Each map should be in A3 size. To be able to create such maps, what libraries or programming languages should I use? I have experience in GeoPandas library but I am not sure if it is enough for this task. Should I try QGIS plugin development or ArcPY? Could you share your experiences and ideas about this problem?


texts:
My clients have an existing ESRI ArcGIS Server with an old Silverlight API based online mapping tool.

As part of a large development program a non ESRI tool will be used to host many forms of data resulting from a large range of automated processes, only some of which will have Lat/Long components.

Is it possible to register an external (non ESRI based) OGC server with ArcGIS Server, so a single restful query to ArcGIS server via say an identify tool will result in both the ArcGIS Server based, and the external OGC based data being queried and then returned?


texts: Display multiple ASCII (xyz) files in QGIS

texts:
I am new to QGIS and I am trying to import more than 700 xyz (it is the same as ASCII?) files into the map so I can get a digital elevation model and then get the average elevation for some small polygons. I can simply import individual xyz file but importing more at once doesn't work for me?

I would like to make something such as a virtual raster catalog but I get:

  Warning 6: gdalbuildvrt does not support positive NS resolution.

My files look like this: (they have no headers)

383000.0  97000.0   83.30
383005.0  97000.0   82.80
383010.0  97000.0   82.70...


texts: Passing ModelBuilder ouput as new variable in ArcPy script?

texts:
I am trying to connect three models I build in ModelBuilder within an ArcPy script.

The first model gets its parameters defined by the user. Whereas the other two models will get some parameters from the output of the first model. I use inline variable substitutions in all of the models which leads to the main problem.

The output of the first model (.shp) is saved as a python variable that I try to pass as a parameter to the second model. When I print this variable, I get the filename with the inline variable substitutions (%Name%_filename.shp).

print(kom_ug)
# C:/Users/bjbr/Desktop/test_bjoern/test/shape\%Name%_ug.shp

This causes the error. The second model isn´t able to find the variable in the directory since the shapefile is stored as file_filename.shp.

Runtime error  Traceback (most recent call last):   File "<string>", line 1, in <module>   File "C:\Users\bjbr\Documents\ArcGIS\Toolbox.tbx", line 26, in KomBufferTool     ExecuteError: Failed to execute. Parameters are not valid. ERROR 000732: Kommune_Endringer_ug: Dataset C:/Users/bjbr/Desktop/test_bjoern/test/shape\%Name%_ug.shp does not exist or is not supported The value does not exist.

When I hardcode the variable name the error disappears, but I would like to work with these inline variable substitutions to comfortably handle a higher volume of data.

Does someone have an idea how to work with inline variable substitutions when calling a model in ArcPy?


texts: WMS GetFeatureInfo and CQL on Multiple Layers

texts:
I am having issues with querying CQL filtered layers with OpenLayers and GeoServer, using the GetFeatureInfo control. It is ignoring the CQL filter and just querying the entire layer (as evidenced by both the results and the request it sends to GeoServer).

I've seen examples where the CQL filter is passed to GeoServer along with the GetFeatureInfo request, but I have several layers visible and each of them has a different CQL filter; so I'm not sure how/if I can handle this by appending my filter to the request URL.

Does anyone have any ideas how I can solve or work around this issue?


texts: Why does not the coordinate of the graph match the dataframe?

texts:
I drawing the plot of Île-de-France with
library(sf)
library(ggplot2)

df <- read.csv("https://raw.githubusercontent.com/leanhdung1994/Stat_Cons/main/station.csv")[1:3]
df2 <- st_as_sf(df, coords = c('V1','V2'), crs = 2154)

df3 <- ggplot() +
  geom_sf(data = df2,
          size = 0.1,
          color = "lawngreen",
          fill = "lawngreen") +
  ggtitle("Stations")

plot(df3)

My dataframe is
Simple feature collection with 25338 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 48.66986 ymin: 2.012258 xmax: 49.05865 ymax: 2.782492
projected CRS:  RGF93 / Lambert-93
First 10 features:
    X                  geometry
1   1 POINT (48.85153 2.637754)
2   2 POINT (48.83918 2.638771)
3   3 POINT (48.85163 2.621466)
4   4  POINT (48.9135 2.482264)
5   5 POINT (48.99439 2.523624)
6   6 POINT (48.88744 2.325725)
7   7 POINT (48.85602 2.624391)
8   8  POINT (48.8531 2.281198)
9   9 POINT (49.00259 2.535003)
10 10 POINT (48.84394 2.627361)

My plot is

Could you explain why the coordinate from the plot is somewhere 5.98°S 1.7°W? I expect it would be somewhat 48.85°N 2.63°E.


texts: ArcMap toolbox script with draw geometries parameter

texts:
I have an arcpy script in toolbox, that need one additional parameter to draw geometry(ies) which will work as a filter. I have done something similar using arcobjects using GPFeatureRecordSetLayerClass as an input parameter and creating temporary layer. If feature set exist, then ArcMap allow to draw geometries in it. Something similar is done if default ArcMap tools "Server Tools>Data Extraction>Extract Data" where feature set "schema" parameter have a dynamic value [toolbox folder path][shapefile] that look something like this

C:\Users\... ...\Desktop10.1\ArcToolbox\My Toolboxes\aoi.shp

When I tried to insert scheme path manually, it tells me that files does not exist. Any ideas to how to replicate feature set behavior from the default "Extract Data" tool?


texts: Exporting cloud masked images from MODIS EVI collection as separate images

texts:
I am performing cloud masking on a collection of modis images.The code below  works on the cloud masking part however when i export the images it brings the error below    .I do not know what is the problem causing the error.

var maskEmptyPixels = function(image) {
      var withObs = image.select('EVI').gt(0);
      return image.updateMask(withObs);
      };
    var maskClouds = function(image) {
      var QA = image.select('SummaryQA');
      var bitMask = 1 << 10;
      return image.updateMask(QA.bitwiseAnd(bitMask).eq(0));
    };
    var table_bounds = function(image){
     return image.clip(table);
    };
    var collection =ee.ImageCollection('MODIS/006/MOD13Q1' )
                   .filterDate('2019-03-01', '2019-07-30')
                    .filterBounds(table).select('EVI','SummaryQA')
                    .map(maskEmptyPixels);
    var evicollection = collection.map(table_bounds);
    var totalObsCount = evicollection
            .select('EVI')
            .count();
    var collectionCloudMasked = evicollection.map(maskClouds);
    var clearObsCount = collectionCloudMasked
            .select('EVI')
            .count()
            .unmask(0);
    Map.setCenter(35.94,-0.37,8);
    Map.addLayer(
        collectionCloudMasked.median(),
        {bands: ['EVI'],
         gain: 0.07,
         gamma: 1.4
        },
        'median of masked collection'
      );
      print(collectionCloudMasked);
    var list = collectionCloudMasked.toList(11);
    for(var i=1;i<11;i++){
      var nam=ee.String('EVI_img')
                 .cat(ee.String(ee.Number(i)))
                 .getInfo();
    }
    Export.image.toDrive({
           image: collectionCloudMasked,
           description: nam,
           scale:250,
           maxPixels:1e13
    });


texts: Cannot acess application created with Python scripting extension of GeoServer?

texts:
I am trying to use the python scripting extension of Geoserver version 2.13.0 to create an application.

I installed the extension by following this tutorial: http://docs.geoserver.org/stable/en/user/community/scripting/installation.html

I created a new script from menu >scripts >add new script

The code I am using :

def app(environ, start_response):
    start_response('200 OK', [('Content-Type', 'text/plain')])
    return ['Hello World!']

My problem is when I try to access the application from the link
ip-address/script/apps/name_application I get the 404 error.

Also when I go to the link ip-address/script/apps I get:

{"scripts": ["name_application"]}

I am following this tutorial:
http://docs.geoserver.org/stable/en/user/community/scripting/hooks.html

Is there anything I am missing here?

A link to Geoserver log file: https://pastebin.com/DfyUJfAc


texts: Added CAD file to ArcMap, annotations that are present in AutoCAD showing up as points in ArcMap

texts:
So I added this cad file to my map using the add data button and everything showed up, except for these texts that belong to the lines. The text itself showed up as points
The annotation layer transferred over no problem but these text values are separate from that annotation layer. The person that made the AutoCAD map only made these details available from what I found in the 'value' column in the attribute editor. From what I could tell, these values did not make it to any of my attribute tables in ArcMap. What would be the best way to migrate that data over?


texts:
I have a TIFF file which I need to convert to dataframe and back to TIFF. This is what I did so far:
from osgeo import gdal
import pandas as pd
import numpy as np

ds = gdal.Open('t2m_20210523.tif')

ar = ds.GetRasterBand(1).ReadAsArray()
flat = ar.flatten()
gt = ds.GetGeoTransform()
res = gt[1]
xmin = gt[0]
ymax = gt[3]
xsize = ds.RasterXSize
ysize = ds.RasterYSize
xstart = xmin +res/2
ystart = ymax - res/2
ds = None
x = np.arange(xstart, xstart+xsize*res, res)
y = np.arange(ystart, ystart-ysize*res, -res)
x = np.tile(x, ysize)
y = np.repeat(y, xsize)

dfn = pd.DataFrame({"x":x[0:len(flat)], "y":y[0:len(flat)], "value":flat})
dfn['value']=dfn['value'].fillna(-999)
dfn=dfn.sort_values(['y','x'])

dfn.to_csv("df.xyz", index = False, header = None, sep = " ")
demn = gdal.Translate("out_file.tif", "df.xyz", outputSRS = "EPSG:4326")
demn = None

This is what I got:

I believe I somehow flipped coordinates, but so far I failed to find combination which works.


texts:
I have scripted tools in ArcPy that gather inputs using "GetParameterAsText", which then show up in the tool GUI like this:

But that is for gathering a single file by navigating to it and selecting it.

I want to script a tool that allows the user to drag multiple layers from the Table of Contents, just like the "Batch Project" tool below.

What is the proper scripting syntax to generate the same input field as "Input Feature Class or Dataset", and also the large field below it where the multiple layers are displayed?  I am using Arc version 10.2.


texts:
I can't make this work. I have a JavaScript program that should upload a KML and it does not work. Uploading the same KML manually in Google Maps works perfectly

<!DOCTYPE html>

    <meta name="viewport" content="initial-scale=1.0, user-scalable=no">
    <meta charset="utf-8">
    <title>GeoRSS Layers</title>


texts: R sf count number of Polygons in Multipolygon

texts:
How can I count the number of polygons that are present in a multipolygon with the R package sf? Example: In the following image that represents a multipolygon I count 3 polygons. I want to be able to store this information as an attribute. Below I added the code that creates the multipolygon.

# Create Multipolygon: https://cran.r-project.org/web/packages/sf/vignettes/sf1.html
library(sf)
p1 <- rbind(c(0,0), c(1,0), c(3,2), c(2,4), c(1,4), c(0,0))
p2 <- rbind(c(1,1), c(1,2), c(2,2), c(1,1))
pol <-st_polygon(list(p1,p2))
p3 <- rbind(c(3,0), c(4,0), c(4,1), c(3,1), c(3,0))
p4 <- rbind(c(3.3,0.3), c(3.8,0.3), c(3.8,0.8), c(3.3,0.8), c(3.3,0.3))[5:1,]
p5 <- rbind(c(3,3), c(4,2), c(4,3), c(3,3))
mpol <- st_multipolygon(list(list(p1,p2), list(p3,p4), list(p5)))

mpol_geom = st_sfc(mpol, crs = 4326)

mpol_attrib = data.frame(
name = "MULTIPOLYGON"
)

mpol_sf = st_sf(mpol_attrib, geometry = mpol_geom)

plot(mpol_sf)


texts:
I am attempting to acquire imagery from the Worldview-2 satellite, and finding limited coverage for many areas within the United States. Sometimes, coverage is temporally deep for a time period (say 15 June 2018 - 31 October 2018) in limited areas, but no coverage for other areas.

Is this a common occurrence?

I am wondering if I am doing something wrong when trying to download the data.


texts:
Let's say that green polygon is under red and part of blue with yellow stripes is also under red polygon. How can I make one polygon from every polygons that intersects my main ( red ) polygon and also remove polygons ( all green and part of blue with yellow stripes ) under red polygon?

I'm using QuantumGIS and PostGIS.

I did this (in QGIS):

dissolve
multipart to singleparts

Please tell me that I'm right ;)


texts: I would like to install QGIS on a cloud platform, like AWS/Amazon cloud for example. Do you think this is a valid solution? For my work I have to upload several satellite images (Sentinel and Landsat) all at a time and perform some processing. I am wondering if I should buy a good PC or if working in the cloud is better.


texts: Redefining the border of a polygon in Postgis

texts:
I have a PostGIS algorithm that writes a table with "pixelated" polygons, that I will call "clusters".

My idea is to smooth the edges of these clusters using concave hulls (ST_ConcaveHull). To do that, I will have to decompose the current edges into their points (using ST_DumpPoints?) and collect them again (ST_Collect?) in order to send them to the concave hull function.
My difficulty is how to produce a query, that in one go can perform these operations, and make sure that the concave hull constructor is called as many times as the clusters, having as argument the collection of points for each cluster.

Something like this:

SELECT ST_Concave_Hull(st_collect((st_dumppoints(geom)).geom),0.8) from exploded group by (st_dumppoints(geom)).geom;

But it's returning this error:

HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

Any ideas? :-(


texts: QGIS export to PDF - rasterize only basemap layers but keep others as vectors

texts:
I am building my own map in QGIS 3.4 in which I created a basemap from ESRI satellite imagery and a ESRI hillshade overlay. On top I have several custom vector layers (country borders, city names etc.)

When exporting the map to PDF, I would like to keep the vector information of labels, borders, lines etcetera to ensure sharpness. The background basemap layers however can (and should) be merged rasterized. Is there a way to only rasterize certain layers to when exporting?


texts: CreateFeatureclass_management doesn't add a layer to the TOC of ArcMap10

texts:
I created a python script for a toolbox and I am wondering about the fact that CreateFeatureclass_management doesn't add the feature-class as layer in the actual dataframe. This does not happen in the environment of a toolbox, but if I create a feature-class in the python window of ArcMap it creates the wanted layer.

I already tried MakeFeatureLayer_management followed by AddLayer, with no result.

What I am doing wrong?


texts: How can I merge thousands of polygons in Oracle Spatial in faster way?

texts:
Previously I was using  sdo_aggr_union(sdoaggrtype(geom, 0.005)) to get boundary polygon by merging geometries, which was taking at least 3 minutes in case of 4000 polygons.

Even if I use merging of polygons by choosing block by block of polygons set, its still taking 30 seconds. (Refer Oracle Complex spatial queries: Examples C-4).

Is there any other approach to achieve such result in minimum time?


texts: Generating vector geospatial PDF from ArcGIS Server services

texts:
I am exploring possibilities to generate vector geospatial PDF from data which is accessible from ArcGIS Server services. I have ArcGIS Server MapService with complex styles. They are very simplified in FeatureService - so that is not an option.

I have already tried to use Export Map operation of MapService and was lucky to get nicely styled vector PDF from it. Of course, I have coordinates of a BoundBox (this info actually is used for querying).

So, now I need to include additional data (some labels, borders and other pre-print graphic) and generate result PDF. It seems like a task for FME Desktop. But I have not found how I can import input PDF from http query. In fact it seems that I can't even open PDF from disk in FME; it is capable only to write it.

What possible solutions do you see for this task?

The key thing is to keep complex styling from ArcGIS Server. I can also access the mxd which was used for creation of the service.

I need solution which is possible to run from scripts, as it will be scheduled task. Maybe there is a solution to combine several PDF files into one in Python? Some advanced Python lib for editing PDFs?


texts: ECW support on ubuntu 14.04 - QGIS 2.8.1

texts:
I've followed instructions for lubuntu 14.04 - QGIS 2.2 ( Answer in Unable to install ECW support on lubuntu 14.04 ), adapting to QGIS 2.8.1. No problem to install packets, but result of gdalinfo --formats | grep -i ECW  is:

ERROR 1: /usr/local/lib/libNCSEcw.so: undefined symbol: _ZTI10CNCSThread
ERROR 1: /usr/local/lib/libNCSEcw.so: undefined symbol: _ZTI10CNCSThread

And QGIS doesn't admit ECW still.
How can I fix it?


texts: ArcGIS JS API: Zooming with mouse wheel/double-click forces pan

texts:
I have a map that I'm initializing off-screen and sliding in using jQuery. For some reason zooming in with the scroll wheel or by double-clicking causes the map to simultaneously pan down (or pan up when zooming out). Using the slider controls works fine. Also if I disable the animation and create it on-screen it zooms normally. I'm guessing it has something to do with the relative position of the map to where it was created, but I don't imagine I have much control over that. Does anyone have any thoughts on this?

Here's how I'm styling the map:

#map {
    bottom: -450px;
    height: 450px;
    position: fixed;
    width: 100%;
}

And this is my animation function:

toggleMap: function() {
    // New bottom value
    var newBottom;

    if (isShowing) {
        newBottom = -450px;
    }
    else {
        newBottom = 0;
    }

    isShowing = !isShowing;
    $map.animate({
        bottom: newBottom
    });
},


texts:
I am able to draw and create wms mapfile for POLYGON and LINE but not for a POINT, I am not sure what maybe the issue, can you assist

LAYER
    NAME    Country
    TYPE POINT
    EXTENT -25 -35 51 37
    CONNECTIONTYPE  POSTGIS
    CONNECTION  'host=localhost dbname=xxx user=xxxx password=xxx port=5432'
    DATA "geom from (SELECT geom, gid FROM education WHERE country = 'xxxx' AND electricity = 'Yes' AND water = 'Yes') as subquery USING unique gid USING srid=4326"
    METADATA
        'wms_title' 'Country'
        'wms_server_version'    '1.1.0'
        'gml_include_items' 'all'
        'wms_include_items' 'all'
    END
    STATUS  OFF
    TRANSPARENCY 100
    PROJECTION
        'init=epsg:4326'
    END

    CLASS
        NAME    'Country'
        STYLE
            SYMBOL "circle"
            SIZE 16
            COLOR 255 0 0
            OUTLINECOLOR 0 0 0
        END
    END
END

Please assist....


texts: Geoserver MBTiles plugins performance issues

texts:
I am using the latest stable 2.7.2 version of GeoServer and I am trying to play with the MBTiles plugin, I have installed also WPS plugin that is needed by the MBTiles one (otherwise I got error when starting). I am using the MBTiles plugin found on community-latest because no one is found under the release folder for 2.7.2 version, however it seems to work, the only problem is that is really bad in performance (30-40s or even more for each request to be displayed). Reading MBTiles under OpenMap (BBN Technologies) or under World Wind Java (I have wrote a plugin for that) is really fast and straight forward, I don't understand why it is so slow. Could be probably something related to different projection system? My MBTiles is in EPSG 3857, but requested display projection is the same, so I suppose that no re-projection is needed....Prior to investigate on the GeoServer code for the MBTiles plugin, I would ask you if someone has experienced the same performance problem.

Taking a thread stack I see that Geoserver seems always executing this piece of code:

at org.sqlite.jdbc3.JDBC3PreparedStatement.executeQuery(JDBC3PreparedStatement.java:68)
at org.apache.commons.dbcp.DelegatingPreparedStatement.executeQuery(DelegatingPreparedStatement.java:96)
at org.apache.commons.dbcp.DelegatingPreparedStatement.executeQuery(DelegatingPreparedStatement.java:96)
at org.geotools.mbtiles.MBTilesFile.closestZoom(MBTilesFile.java:519)
at org.geotools.mbtiles.mosaic.MBTilesReader.read(MBTilesReader.java:152)
at org.geotools.coverage.grid.io.AbstractGridCoverage2DReader.read(AbstractGridCoverage2DReader.java:249)
at org.geoserver.catalog.SingleGridCoverage2DReader.read(SingleGridCoverage2DReader.java:147)
at org.geoserver.catalog.CoverageDimensionCustomizerReader.read(CoverageDimensionCustomizerReader.java:227)
at org.geoserver.catalog.CoverageDimensionCustomizerReader.read(CoverageDimensionCustomizerReader.java:217)


texts:
I want to run the spatial join tool to join hundreds of shapefiles from two different folders.

For example one folder contains point  shapefiles with names d1.1,d1.2,d1.3  . . . .

the other folder contain shapefiles with names q1.1,q1.2,q1.3  . . . .

i want to join d1.1 to q1.1, d1.2 to q1.2  and so on, using python but getting the error:

Traceback (most recent call last):

  File "C:\spatial join.py", line 28, in <module>

    arcpy.SpatialJoin_analysis(tarf, joinf, output, "JOIN_ONE_TO_ONE",

"KEEP_ALL", "POINTID \"POINTID\" true true false 6 Long 0 6 ,First,#,fc[:-
4],POINTID,-1,-1;AOD \"AOD\" true true false 10 Long 0 10 ,First,#,fc[:-4],GRID_CODE,-1,-1;QF \"QF\" true true false 10 Long 0 10 ,First,#,fs[:-4],GRID_CODE,-1,-1", "INTERSECT", "", "")

  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\analysis.py", line 471, in SpatialJoin
    raise e

ExecuteError: Failed to execute. Parameters are not valid.

ERROR 000732: Target Features: Dataset da2016001.0510 does not exist or is not supported

ERROR 000732: Join Features: Dataset qfa2016001.0510 does not exist or is not supported

Failed to execute (SpatialJoin).

Here is my code

import arcpy

arcpy.env.workspace = r"C:\pointsD"

fcsD = arcpy.ListFeatureClasses()

arcpy.env.workspace = r"C:\pointsQF"

fcsQF = arcpy.ListFeatureClasses()

for fc in fcsD:

    for fs in fcsQF:

        if fs[3:]== fc[2:]:

            tarf = fc[:-4]

            joinf = fs[:-4]

            output = r"C:\DWithQF" + fc[1:-4]

            # Process: Spatial Join

            arcpy.SpatialJoin_analysis(tarf, joinf, output, "JOIN_ONE_TO_ONE",
"KEEP_ALL", "POINTID \"POINTID\" true true false 6 Long 0 6 ,First,#,fc[:-4],POINTID,-1,-1;AOD \"AOD\" true true false 10 Long 0 10 ,First,#,fc[:-4],GRID_CODE,-1,-1;QF \"QF\" true true false 10 Long 0 10 ,First,#,fs[:-4],GRID_CODE,-1,-1", "INTERSECT", "", "")

print "Done"


texts: Intersection and add selection in pyqgis

texts:
I'm trying to add the intersected geometries of 2 layers in a pyqgis code.

I have a grid made up by cells (so a polygon) and a line layer. I want to select the features of the grid layer where they intersect the line features (so basically the Select by Location tool used within a python script).

The code (that works) is the following:

# get the layers from the interface
grid = iface.activeLayer()
line = iface.activeLayer()

for i in grid.getFeatures():
    for l in line.getFeatures():
        if i.geometry().intersects(l.geometry()):
            print i

It returns (correctly) 17 features.

What I'm not able to to is to add these QgsFeatures as a selection to the grid layer.


texts: qgis2web plugin overriding system locale when publishing (comma as decimal separator)

texts:
I'm using qgis2web to export maps and even though in desktop client my settings are all right, when I publish (with both Leaflet or OpenLayers) the decimal separator changes from comma to point. Is there a way to fix it?

As I said, this happens only when exporting, and my concern is exclusively about how values are shown in popup fields / legend. I've considered converting values to string to overcome this, but exported files are already too big and didn't want to add any new fields (some fields need to remain as real for categorization).

I'm using QGIS 3.8.3 and qgis2web 3.12.1.


texts:
I am using MapProxy to use my company's tiled map service (TMS) as a web map service (WMS) for consumption by desktop GIS clients. One of the great advantages is being able to add additional projections to the configuration YAML and for them to render on the fly; I have been successful in adding support for Antarctic and Arctic Stereographic Polar projections.

Really the only issue left is that MapProxy does not play nice with the antimeridian. When an image is requested that consists of tiles either side of the antimeridian, regardless of projection, rather than requesting the tiles that abut the antimeridian and stitching them together, it seems to request the entire map in a "reacharound" fashion, and then squeeze all of that information into a small part of the map surrounding the antimeridian. A picture can probably explain better:

Ignore the fact that New Zealand is upside-down, it's just a polar projection. If I zoom in, the region of "smear" becomes smaller in terms of geographical spread, but stays roughly the same size on my screen. The opposite applies if  if I zoom out; indeed at one zoom level out from the image, half of New Zealand is affected by the smear. The same applies in the Northern Hemispshere at the Bering Strait and indeed right along the antimeridian in every projection I have tried.

I have tried changing my configuration YAML file in the following ways, not really expecting them to solve the problem but trying it anyway:

Changing the bounding box of the grids from (for example) [-180,45,180,89] to [-179,45,179,89] (I'd rather miss a small range than have a larger range overruled by the image smear). I tried this with the non-default projections I was using, as well as overriding the base grid GLOBAL_MERCATOR bounding box property.
Changing the meta_buffer and meta_size to 0, in case the buffer was somehow causing this.

Has anyone else encountered this issue and has a solution or a workaround?

EDIT: I am using MapProxy v1.8.0


texts: Deleting nearly identical points in one dimention

texts:
I have extracted values from a raster in table form, done some calculations and plotted these points as XY data, in which I hope to use to create a line. However due to the original layout of the raster, I occasionally get a chatter effect. I would like to be able to delete the points where ONLY the X-coord values is within 0.7m of another x-cord value.

I have tried the Integrate tool, however it seems to take into account the Y coordinate as well. I have tried to smooth the line without luck as well.

I have no idea how to approach this problem.

UPDATE
After running the code below in on the line feature class I got the following result, it basically skipped over all the duplicate points, which serves the purpose of what I needed.  I don't know why the last point was dropped. Old line is in black, the new line is in red.

A more elegant approach would create and use a midpoint between each pair of the close points.


texts: featureInfo limited to first match after qgis2web WMS export

texts:
Here's my problem:

I created a WMS using PostGIS and GeoServer and integrated it into QGis. So far everything worked fine. Then I wanted to display the WMS as a webpage and exported it using qgis2web.

As you can see in the picture (right the WMS in QGis, left the WMS as a webpage, information refers to the same polygon) not all information assigned to the polygon is shown in the popup of the webpage. Only the first line of information is displayed.For this particular Polygon about 10 rows of Information are missing.

EDIT:
After Ian Turton helped me, it turned out that the feautureInfo of the qgis2web export of the WMS is limited to the first match.

But my question remains the same:
Is this a problem caused by the export in the WMS, or is it a problem I have in the web page because feautureInfo doen´t read all the information?

Thanks


texts: Server side reporting options for Silverlight Viewer

texts:
I am looking for suggestions from user community regarding some options for complex reporting. I would like to allow users to browse for data and then based on their selection they can generate a preformatted report and download the reprot as PDF. I want to have reporting logic on server so that we can accommodate any future technology changes and be able to use existing development

I have looked at number of third party controls like Infragistics/Comopnent One/Pertpetuum/ Dev Express/Telrik/Active Reports. I am not sure Crystal reports are an option. Most of the third party controls that I have seen are based on client side.

Appreciate any feedback.

Thanks
Jay


texts: saving a SVG marker to a layer for sharing

texts: I have a layer which consists of point data. I have switched the symbology of the layer to a svg image that I loaded in. To maintain the symbology I went to LayerProperties>style and then i saved the style. However, when I open the layer on a different computer (one with out the svg saved in) and apply the layer style nothing come up. How can I save a picture marker to a layer file and still be able to send it via email etc. and have it show up on another computer.


texts:
I can see that for the Aqua instruments (AIRS, MODIS, CERES etc), the Level 2 data are binned into Level 3 data providing average values for (say) a 1° x 1° cells grid.
Let's say that I would like to bin the data into a different size grid (e.g. 0.5° x 0.5°). I know that I could just write some code that takes Level 2 data as input and provides the binned data as output, but I wonder if there are tools that do this job.
I have found the pymodis python library and the SeaDAS l2bin program which might be useful but I'm not sure I'm on the right track.
Any suggestion ?


texts: Creating extent indicator in Data View of ArcMap?

texts: I know how to create an extent indicator in Layout View, but I have no idea how to create one for Data View. I have 25 bookmarks created for Data View and what I need to do is to display all of them in Data View. I work in ArcMap 10.3.


texts: Changing Malvinas/Falkands (Argentina, South America) label in CartoDB?

texts:
I'd like to know if there any way to change or hide the label of "Falklands Islands" in south Atlantic Ocean, South America, Argentina, in CartoDB.

I know CartoDB works with Open Street View.


texts: ArcGIS Runtime SDK for Java - WebMap display issue (spatial reference construction error)

texts:
I'm trying to display a map which I created and uploaded on ArcGIS online as a map package:

map = new JMap();
Portal portal = new Portal("http://www.arcgis.com", null);
PortalItem item = null;
item = PortalItem.fetchItem(portal, "9464b884blahblahblah");
WebMap wm = WebMap.newInstance(item);
map.loadWebMap(wm);

I tried the code with a map ID from the samples provided in the docs and it works fine.
With the ID of the map I uploaded, however, I get this error:

Java version : 1.8.0_20 (Oracle Corporation) amd64
Rendering engine : DirectX
Unable to construct spatial reference from GCS_WGS_1984
org.codehaus.jackson.JsonParseException: Unexpected character ('z' (code 122)): expected a valid value (number, String, array, object, 'true', 'false' or 'null')
 at [Source: org.apache.http.conn.EofSensorInputStream@48718cc6; line: 1, column: 3]
    at org.codehaus.jackson.JsonParser._constructError(JsonParser.java:1432)
    at org.codehaus.jackson.impl.JsonParserMinimalBase._reportError(JsonParserMinimalBase.java:521)
    at org.codehaus.jackson.impl.JsonParserMinimalBase._reportUnexpectedChar(JsonParserMinimalBase.java:442)
    at org.codehaus.jackson.impl.Utf8StreamParser._handleUnexpectedValue(Utf8StreamParser.java:2090)
    at org.codehaus.jackson.impl.Utf8StreamParser._nextTokenNotInObject(Utf8StreamParser.java:606)
    at org.codehaus.jackson.impl.Utf8StreamParser.nextToken(Utf8StreamParser.java:492)
    at com.esri.core.portal.WebMap.a(Unknown Source)
    at com.esri.core.portal.WebMap.<init>(Unknown Source)
    at com.esri.core.portal.WebMap.<init>(Unknown Source)
    at com.esri.core.portal.WebMap.newInstance(Unknown Source)
    at Tsunami.<init>(Tsunami.java:80)
    at Tsunami$2.run(Tsunami.java:105)
    at java.awt.event.InvocationEvent.dispatch(InvocationEvent.java:311)
    at java.awt.EventQueue.dispatchEventImpl(EventQueue.java:744)
    at java.awt.EventQueue.access$400(EventQueue.java:97)
    at java.awt.EventQueue$3.run(EventQueue.java:697)
    at java.awt.EventQueue$3.run(EventQueue.java:691)
    at java.security.AccessController.doPrivileged(Native Method)
    at java.security.ProtectionDomain$1.doIntersectionPrivilege(ProtectionDomain.java:75)
    at java.awt.EventQueue.dispatchEvent(EventQueue.java:714)
    at java.awt.EventDispatchThread.pumpOneEventForFilters(EventDispatchThread.java:201)
    at java.awt.EventDispatchThread.pumpEventsForFilter(EventDispatchThread.java:116)
    at java.awt.EventDispatchThread.pumpEventsForHierarchy(EventDispatchThread.java:105)
    at java.awt.EventDispatchThread.pumpEvents(EventDispatchThread.java:101)
    at java.awt.EventDispatchThread.pumpEvents(EventDispatchThread.java:93)
    at java.awt.EventDispatchThread.run(EventDispatchThread.java:82)

The "Unable to construct spatial reference ..." error message pops up after 1-2 seconds after starting the program, where as the JsonParseException takes about 10 seconds to appear and is triggered by the following line:

    WebMap wm = WebMap.newInstance(item);

I checked the layers I created (in properties, in the coordinate system tab) they do reside under a "GCS_WGS_1984" coordinate system.

Is the JsonParseException somehow linked to the Unable to construct... error message?

Can anybody tell me what am I doing wrong?


texts:
I am looking for the most accurate way of calculating the view angle (Azimuth, Elevation) from a reference point (my location) A to another point B. I have the Lat, Lon, and Altitude for each point.

Currently I am converting to local ENU coordinates using the Wikipedia formulas for ECEF to ENU and using that to calculate the angles.

Az = atan2(E,N)
El = atan(U / sqrt(E^2 + N^2)

Is there a more accurate way of calculating these Angles?

Maybe a different coordinate system I should be projecting to?

Is this the correct type of Angle, should a rhumb line azimuth or something else be used instead?

I have also looked at using Vincenty's formula for calculating the azimuth but it does not take into account the altitude of each point. A change in Altitude using ENU seems to slightly change the Azimuth which seems correct(?). When comparing Vincenty Azimuth with the ENU Azimuth on the surface (Altitude 0), the answers are extremely close but start to differ if the points are not on the surface.