texts:
I'm trying find libraries for R that would find the shortest route between two locations, given a SpatialLinesDataFrame representing the road network. All I could find was stplanr, which delegates the routing to external services such as CycleStreets.net (which does not cover the area I'm interested in) or GraphHopper (which is expensive).

What am I missing? Are there open routing packages for R available on GitHub or CRAN?

(My preference on an R package is due to the routing solution being a part of a larger application, and a native lib would be easier to integrate.)


texts: Trying to add lat long points to a map, not getting an ok button

texts:
I add "Data Source Manager" then "Delimited Text" and i select the x y and coordinate system but there is no OK button?

just add and close add just restarts the process. how do i get the ok button?


texts:
I am trying to convert a .bin file of sea ice concentration to a GeoTIFF by following these directions from NSIDC: https://nsidc.org/support/how/how-do-i-convert-nsidc-0051-sea-ice-concentration-data-binary-geotiff

What I am running in terminal looks like this:

MacBook-Air-2:converting rbex$ gdal_translate -of GTiff -a_srs '+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0+y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs' -a_nodata 255 -A_ullr -3850000.0 5850000.0 3750000.0 -5350000.0 nt_20160101_f17_v1.1_n.bin nt_20160101_f17_v1.1_n.tif

However, when I run the provided code, it returns:

ERROR 4: `nt_20160101_f17_v1.1_n.bin' not recognized as a supported file format.

I suspect this is because the binary file, nt_20160101_f17_v1.1_n.bin, has a zipper icon on it and so is a zipped file of some sort? I can not find any documentation anywhere on how to unzip the .bin file, or any other kind of trouble shooting to do.

Update:

I did write a header, formatted like this:

ENVI
description = {nt_20160101_f17_v1.1_n.bin}
samples = 304
lines   = 448
bands   = 1
header offset = 300
file type = ENVI Standard
data type = 1
interleave = bsq
byte order = 0
map info = {Polar Stereographic, 1, 1, -3850000, 5850000, 25000, 25000} projection info = {31, 6378273, 6356889.449, 70, -45, 0, 0, Polar Stereographic} coordinate system string = {PROJCS["Stereographic_North_Pole",GEOGCS["GCS_unnamed ellipse",DATUM["D_unknown",SPHEROID["Unknown",6378273,298.279411123064]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Stereographic_North_Pole"],PARAMETER["standard_parallel_1",70],PARAMETER["central_meridian",-45],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]}
band names = {Band 1}

The naming for the header is nt_20160101_f17_v1.1_n.bin.hdr , it is saved in the same folder as the .bin file, and that file is the directory in my terminal.

I tried using Stuffit Expander on nt_20160101_f17_v1.1_n.bin, but the expander says "unable to determine the file type." When I try unpacking using my mac, it creates a file that ends in .bin.cpgz

My data is coming from: https://daacdata.apps.nsidc.org/pub/DATASETS/nsidc0051_gsfc_nasateam_seaice/final-gsfc/north/daily/2015/
(this requires an earth data account to access)


texts: Layer of roads segments with attribute "STEL" by rsd.cz NDIC

texts:
Is it possible download layer with attribute "STEL"?
Attribute means segments of roads in Czech republic created by "ŘSD", more information at this link.
"STEL" attribute is using for visualization NDIC data, more informations about NDIC are at this link.
Documentation for NDIC (english language):
link
Documentation for NDIC (czech language):
link


texts: Using convex hull for interpolation of individual plots or lump plots together in ArcMap?

texts:
This question is inspired by Using kriging to extrapolate values outside of sampling polygon in ArcMap?

I'm working on a project trying to understand the lasting effects of historic scientific activities on delicate soil communities. Essentially, the hypothesis is that there will be fewer soil organisms the closer one gets to the center of where human activity took place (due to trampling and other disturbances).

I've got soil data that were collected at known locations of historic camps where researchers lived and worked in the field. After locating the footprint of where each structure in the camp used to be (e.g. hut, shed, etc.), a sample was taken at the center of the footprint and at three points along three different transects radiating out from the plot center.

Should I use the convex hull of each individual plot to interpolate values, or can I do a convex hull and interpolation of the entire camp (based on all sample points)?

I have included simple visual examples of these two options below. The scale is approximately 150m from one end of the camp to the other.

A camp with three plots and their respective convex hulls

Plot-scale samples and interpolation (IDW) (only one plot shown)

Camp-scale samples and convex hull

Camp-scale samples and test interpolation (IDW)

Note: the interpolated surfaces are not scaled to each other in these visuals.


texts: Applying same color ramp to multiple layers in ArcMap (in one step)?

texts: I attached the image where you can see the layers in my tabnle of contents. The first one has the wanted color ramp, which I want to apply to the rest down there. There must be a way to avoid doing it individually on each layer ?


texts:
Can anyone tell me how to import OSM data in mongoDB step by step ?

This is the link to the data http://download.geofabrik.de/europe/germany/baden-wuerttemberg/karlsruhe-regbez.html

https://gis.stackexchange.com/users/12320/ylka

this is the error i am getting in the 1 step:

  Unable to open datasource with the following drivers ... [huge list]


texts: Checking whether wms layer has any features?

texts:
I've been searching the net for a couple of days to find a way, how to check if a wms layer has no features...

I'm using openLayers/geoserver/postgis.
I'm using a cql_filter to search a wms layer. If the filter returns some features (points or lines or polygons) then i redraw the layer and they succesfully render.

But, if no features are rendered, I want to alert a message using javascript. The message simply gives some info, like "no matches found, please try again".

I guess I have to check if there are no features, and then alert the message. Right? Or is there another way to do this?

I was thinking, using a simple getFeatureInfo, with no feature.count parameter, after I redraw my layer. If the text that returns is empty, then there are no features rendered (= no matches found). I came up with this idea, playing with the demo requests on Geoserver, using the getFeatureInfo. The syntax they use there, returns all the features...So, I will do the same thing and if no features returned I get what I want...
But how to check for an empty text? Or this idea is completly wrong?


texts:
I've got a mysql table with points and polygons spatial geometries.

I can see that MySQL 5.5 doesn't have the ST_AsGeoJSON() function to extract to GeoJSON.

I'm looking for the simplest way to load these geometries to Google Maps API.

Before starting to read data in PHP and put them in an array structure and send it in GeoJSON format, I wonder if there is a simplest way to do it.

Any suggestion?


texts: Can extracting raster values work in polygons' CRS without raster resampling?

texts:
I am trying to extract raster values using polygon layer with different CRS. Raster is in WGS84 (EPSG:4326) and the polygons in UTM (EPSG:23031).

1. If I want to work in the projection of the polygons, do the cells of the raster need to be resampled?

I.e. is converting to different pixel grid, and thus information loss, necessary? Or can the GIS just reproject the edges of the raster cells and handle it without resampling?

2. Is there any general/prefered practice that the polygons are converted to the projection of the raster?

E.g. to avoid raster resampling - if answer to question 1 is yes.

In other words, can we work in polygons' CRS without the need for raster resampling?

PS: my use case - I am using raster library in R to do the job. R decides to automatically transform polygons:

> a <- extract(raster, polygons, weights = TRUE)
Warning message:
In .local(x, y, ...) :
  Transforming SpatialPolygons to the CRS of the Raster

However, this means that the already projected polygons (which were actually created as a buffer in UTM projected plane) need to be converted back to WGS84 implicit latlong "projection"... I would prefer to work in the plane defined by UTM, because the projection is much less distorted for my purposes, e.g. weighting by areas would work with minimal distortion etc...


texts: How does a spatial database display maps through WMS?

texts: If a spatial database is configured to display maps through WMS, is it always configured to display through TMS and WMTS also? Or is this something that must be configured on the server side separately?


texts: using Map Service or Feature Service for Feature Layer

texts:
Both, Dynamic Map Services and Feature Services allow for querying data (through FeatureLayer) and bringing FeatureSets to the browser.

Which of the above would yield better performance at querying attributes and geometries of large datasets?


texts: Using clip tool in qgis - contours to exact map area

texts:
I am using QGIS. I am trying to clip contours and water vector layers to the exact size of my map (polygon layer). Sometimes it works; sometimes the program crashes before completion (at the same point each time); and sometimes it runs but doesn't clip the exact area. How can I prevent this from happening?
Thanks


texts: Turning off labels that have null values in ArcGIS Desktop?

texts:
How do I turn off labels that have null values?

I am using a leader with this label but I want it turned off if the values are null


texts: How to create a 1:n join and concatenate text attributes?

texts:
I have 2 polygon layers A and B.
I would like to join the attributes from B into A for all the polygons from B that intersect A and for text attributes, i would like them concatenate.
But apparently "Join attributes by location" only works for numerical attributes.


texts: Why does the same UTM northing give different values when converted to latitude?

texts:
I'm converting two coordinates from UTM to latitude/longitude. The coordinates are:
Easting: 657787
Northing: 7795664

and
Easting: 658144
Northing: 7795664

Both from zone 23S (EPSG:31998) or 23K (following this reference). Notice that both have the same northing.
The latitudes given for both are slightly different:
-19.929034598845494 and
-19.92900563505912

I have tried this using R's spTransform function, and also this website. Why aren't they giving me exactly the same latitude, since the northing is the same?


texts: Assigning new values to polygon based on overlapping area with other polygons

texts:
I have two square polygons df1 and df2. Both polygons are assigned a value of 10 (df1) and 20 (df2) in the "content" column. Let's assume, I have a circle with the midpoint 1.5,1.5 and a radius of 1.

How can I assign to this circle a new value "content", based on the overlap with the other polygons?

For example, one third of the circle overlaps with df1 and one third of the circle overlaps with df2. In this case, the circle would then obtain 10 * 1/3 + 20 * 1/3 as a value for "content".

from shapely.geometry import Polygon
import geopandas

polys1 = geopandas.GeoSeries([Polygon([(0,0), (2,0), (2,2), (0,2)])])

polys2 = geopandas.GeoSeries([Polygon([(1,1), (3,1), (3,3), (1,3)])])

df1 = geopandas.GeoDataFrame({'geometry': polys1, 'df1':[1], 'content':[10]})

df2 = geopandas.GeoDataFrame({'geometry': polys2, 'df2':[1], 'content':[20]})

ax = df1.plot(color='red');
df2.plot(ax=ax, color='green', alpha=0.5);


texts:
I'm using geopandas in a Jupyter notebook, working on OSX. I have imported a couple of shapefiles and transformed their CRSs, without any problems. The shapefiles aren't huge - one of the raw .shp files is 30MB, the other is less than 1MB.

I want to use geopandas to calculate the area of the overlap between them, but when I use overlay to start doing this:

df_sussex = gpd.overlay(df_sussex, df_sac, how='intersection')

the Jupyter notebook just hangs. It's been running for over an hour. There are no errors in the console, and the Activity monitor shows it's taking up 98% of CPU. I'm working on a machine with 16GB of RAM - surely it shouldn't take this long? In PostGIS on the same machine it just takes a few seconds.

Is there any way I can debug what's going on?

Both of the dataframes are in EPSG:4326, so I wonder if it could possibly be related to this geopandas issue.


texts: can anyone explain the difference between .qgz and .qgs project files in QGIS or know a good documentation?


texts: How to add XY point and update attributes with ArcPy 10.0

texts:
I would like to create a script tool that takes 3 parameters in ArcMap:

Google maps coordinate (e.g. -33.876846, 151.215117)
text string
text string

The coordinates would define the XY geometry, and the two strings would be the attributes for the new point ("Report" and "Hyperlink")
I am limited to ArcMap 10.0, so I can't use the .da module. I think I have the geometry figured out, but not the attributes.

import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]

#for use in ArcGIS script - gets input coordinates
coordsIn = arcpy.GetParameterAsText(0)

#get Report name
Report = arcpy.GetParameterAsText(1)

#get Report Hyperlink
Hyperlink = arcpy.GetParameterAsText(2)

#delimits coordinates by comma and places them in list "coordsplit"
coordsplit = coordsIn.split(",");

#defines coordinates
X = float(coordsplit[0])
Y = float(coordsplit[1])

#assign feature class to edit
featureClass = "D:\\data\\testing\\testing.gdb\\ReportsGDA"

#create insert cursor
rowInserter = arcpy.InsertCursor(featureClass)

#create update cursor
rowUpdater = arcpy.UpdateCursor(featureClass)

#assign coordinates to point object
pointGeometry = arcpy.Point(X,Y)

# Use the insert cursor to put the point object
# in the feature class
newPoint = rowInserter.newRow()
newPoint.Shape = pointGeometry
rowInserter.insertRow(newPoint)

#update attributes
for feature in rowUpdater:
    if feature.Report == None:
        feature.Report = Report
        feature.Hyperlink = Hyperlink


texts: Importing Shapefile geometry into SQL Server 2016?

texts: What is the most preferred way to import shapefile geometry into MS SQL Server 2016?


texts: Any advantage to use elongated focal window for raster Pearson's Correlation?

texts:
I am performing Pearson's Correlation with 2 rasters (same res, size, extent), and will use a focal window, but is there any advantage or robustness to using a focal window that better reflects the raster extent and shape? e.g.

require(raster)
# make 2 dummy rasters and stack
r1 <- raster(xmn=0,xmx=30,ymn=0,ymx=70,res=1)
r1[] <- sample(1:3,30*70,replace = T)
r2 <- raster(xmn=0,xmx=30,ymn=0,ymx=70,res=1)
r2[] <- sample(1:3,30*70,replace = T)
s <- stack(r1,r2)

# use raster:corLocal to perform pearson test, 3x3 focal window
r.res1 <- corLocal(s[[1]],s[[2]],ngb=5,method="pearson",test=T)
# use raster:corLocal to perform pearson test, 7x3 focal window to reflect dimensions of input rasters
r.res2 <- corLocal(s[[1]],s[[2]],ngb=c(7,3),method="pearson",test=T)

# stack results and plot
res.st <- stack(r.res1, r.res2)
plot(res.st)

Also, does anyone know why Pearson Correlation in spatialEco divides the correlation by n-1? I thought it did not effect r (link here) but it clearly does if you use two identical dummy rasters that should be r=1 all over


texts:
I've got PostGIS installed and would like to create a web application to do some simple edits:

Draw manually (over a Basemap) points, lines and polylines, using snapping
points (points, end of lines, intersections, center of circles)
Delete points, lines and polylines
Rotate a line or group of lines around a point
Move a line or group of lines with a base point
Move polylines and points (manually and to a point given by the coordinates)
Edit position of vertices of polylines and lines (manually and to a point given by the coordinates)
Create buffers around points, lines and polylines

Which is the best way to do this via web? Is there any already-built JavaScript library?


texts: Converting Lambert Conformal Conic (LCC) to UTM

texts:
I have two numpy meshgrids (X, Y both 2D) in Lambert Conformal conic with the following information:

I'd like to directly (or indirectly if not possible) project them to UTM. I have both pyproj and gdal installed, but wasn't able to figure out how to use them. I found solutions for converting from LCC to WGS84 here and here, but couldn't use either of them since I don't have all the information needed.

Even if somehow I can convert from LCC to WGS84, then I can use Python UTM package and utm.from_latlon to get UTM out of it.

I found the following line of code from here:

awips221 = Proj(proj='lcc', R=6371200, lat_1=50, lat_2=50,lon_0=-107, ellps='clrk66')

And then I can use transform from pyproj. In that case, is lat_1 and lat_2 both 25.0 in my case? What would be lon_0 and ellps?

Using gdalinfo we can see lCC and WGS84 projections for the center and four corners of the grid in the picture below:


texts: Finding the nearest line to a point in QGIS,but Not Work

texts:
I saw Finding the nearest line to a point in QGIS?

I wanted to create a program to get the line closest to the selected point.

Here is the code.

        dPos = mouseEvent.pos()
        qpoint = self.get_canvas_point(dPos)
        layer = self.iface.activeLayer()
        features = layer.getFeatures()

        SpInd = QgsSpatialIndex(features)
        nearestSpIndids = SpInd.nearestNeighbor(qpoint,2)

        for nearestid in nearestSpIndids:
            gettingfeatures = layer.getFeatures(QgsFeatureRequest(nearestid))

            for feat in gettingfeatures:
                #feature_check
                prop1 = feat['PEN_Color']
                prop2 = feat['ID']
                root = QTreeWidgetItem()
                root.setText(0,str(prop1) +':' + str(prop2))

                self.dlg.treeWidget.addTopLevelItem(root)

However, even if I select any places, "Specific" lines will be selected.

For example : Regardless of where on the map you choose, you can only get information on the red line.

Is there anything strange about code?


texts: Converting an unprojected coordinate system to a projected one in QGIS 3.10 without ruining data

texts:
I am new to GIS and I think I may have posted on the wrong website:
https://stackoverflow.com/questions/62595452/how-to-convert-an-uprojected-coordinate-system-to-a-projected-one-in-qgis-3-10-w
I created a map in QGIS 3.10 using a non-projected CRS (WGS84 EPSG:4326). I want to display it in 3D Map View, but QGIS tells me that 3Dview does not support uprojected CRS, and I have to switch to a projected CRS. But when I do that, distances get distorted.
I am writing a novel and this is the map of a fictional, island-like location (a microcosm, a metaphor for earth). It is about 5kms x 10kms, currently placed exactly on the South Pole (I know, probably the root of my trouble).
I guess that's where I made a huge blunder. I did not know that QGIS does not support unprojected CRS for 3Dview. If I knew, I would have created everything using a projected CRS and would have placed my imaginary map somewhere less symbolic and more compatible with projected systems. Now, I have created so much data that it is impossible to start anew (the picture shows only a small portion of that data).
I switched to a projected CRS for both the layers and the project, but (obviously) the distances shrank almost 100,000 times. Using some systems (E.g. WGS84/World Mercator EPSG:3395) did not distort the map. The proportions remained fine, but distances got distorted and instead of 5km, I got 5cm! (That's not a microcosm for humans, it's for ants). Some (E.g. ESRI:102037 - South_Pole_Orthographic) distorted both.
I want to convert my unprojected CRS to a projected one, without distorting proportions or distances. Since this is a fictional place, the location of the map, real coordinates, and the CRS are NOT important. As long as the shape and distances remain intact, any location, CRS, etc. will do.
How do I do the conversion with the least hassle?


texts: simple plotting of (latitude, longitude) points on a global map

texts:
I have a dataframe like the following:

country <- c('US', 'CA','TW')
city <- c('Atlanta','Halifax','Taipei')
lat <- c(33.8004, 44.6500, 25.0392 )
lon <- c( -84.3865,-63.6000,121.5250 )

data <- data_frame(country,city, lat, lon)

> data
# A tibble: 3 × 4
  country    city     lat      lon
    <chr>   <chr>   <dbl>    <dbl>
1      US Atlanta 33.8004 -84.3865
2      CA Halifax 44.6500 -63.6000
3      TW  Taipei 25.0392 121.5250

Think of this data as each observations is a mobile call from a city around the world (there could be multiple call from the same city of course).

There are two things I would like to do with this data.

Simplest map. I would simply like to (nicely) plot these calls as a many dots on a global map. So dots may of course be overlaid (with some alpha blending)
Nicer map. First, count how many observations per cities (easy to do with dplyr), and then create a map where polygons are cities. That way, each city has a color according to the number of calls in that city. Ideally I would like to get something like

but for cities on a global map.

Problem is: I dont know where to get the right maps for this and to complicate the matter I can only download the maps manually (not via get_map for instance)

Any ideas how to get these two simple maps?


texts:
Reading the documentation...it looks like I should have this right, but I don't.

The example is:

  Point(0,0).distance(Point(1,1))

  1.4142135623730951

Based on this example I wrote this block (..note: update, 'capital' starts as a shapely obj, and city_pt is a dict...full of multipoints it appears, looking like this).

[{'geometry': {'type': 'MultiPoint',
'coordinates': [(4942585.391221348,  #3940520.723517349)]},
'type': 'Feature', 'id': '17', 'properties':
OrderedDict([(u'Status', u'Bad')])}]

    # capital = shapely obj
    capital_pt = capital.coords
    # city_point is a dict ... multipoints
    for city_point in filtered_all:
        city_items = shape(city_point['geometry'])
        # capital_pt = POINT (13531245.47570414 2886003.268927813)
        # city_items[0] = POINT (4942585.391221348 3940520.723517349)

        # measure distance

        # calculate the distance between capital and cities
        distance_between_pts = capital_pt.distance(city_point)

and got this:

    AttributeError: 'CoordinateSequence' object has no attribute 'distance'

trying the om_henners answer I get

    AttributeError: 'dict' object has no attribute 'coords'

Just in terminal I tried this:

      x = POINT (13531245.47570414 2886003.268927813).distance(POINT (4942585.391221348 3940520.723517349))

      File "", line 1
          x = POINT (13531245.47570414 2886003.268927813).distance(POINT (4942585.391221348 3940520.723517349))
                                                       ^
      SyntaxError: invalid syntax


texts: Problem with Grid Statistics For Polygons

texts:
So honestly I have no idea what I'm doing, I'm a total beginner with GIS and I'm just following directions for an assignment given by my teacher.

I am trying to use Grid Statistics For Polygons -tool to do create statistics(?) of a raster layer

Maybe this image will give a better idea of what I'm doing:

And the log:

processing.runalg(
    "saga:gridstatisticsforpolygons",
    "C:/Users/Name/Pictures/オマタ/H17/Viheralueet.tif",
    "C:/Users/Name/Pictures/オマタ/H17/Bufferi50.shp",
    False, False, False, False, True, False, False, False, 0,
    "C:/Users/Name/Pictures/Viherbufferi.shp")

I can't open the file the tool creates once I've clicked 'Run'. I'm not sure, but is the tool supposed to create a .shp -formatted file? Because the only file the tool creates is Viherbufferi.prj and I'm not able to do anything with it. According to the directions I'm supposed to be able to open the file as a new layer.

I hope someone understood my explanation!


texts: how to digitize areas of a raster and snap to the raster cells?

texts:
I need to digitize the extent of an ecw raster image where the RGB value is not black.(see attatchment)ie. not the envelope, just the area where the image is shown without the border areas. The black areas are not always 0,0,0 which is the main problem. Is there any way to snap to the cell edges?

Available Software:ArcInfo, FME, GDAL


texts: Getting arcpy.da.walk to continue walking after encountering corrupt/problem file

texts:
I am relatively new to python. I am trying to iterate through all DWG polylines in a directory with multiple (hundreds) of subdirectories. I would like the code to skip any DWG files that it cannot do something to and continue iterating.

The following code works fine, until it comes across any DWG files that are possibly corrupted. Python stalls and I cannot get it to throw an exception.

I have isolated the problem to the arcpy.da.walk. It appears that the walk function can open the corrupted DWG to find the polyline, but it cannot get out of the DWG to continue the walk.

I have tried multiple exception types. I am not sure how the onerror function call works in arcpy.da.walk,or even if this is an avenue of exploration (I cannot find any examples where onerror is anything other than "NONE").

import arcpy
import os
import sys
from arcpy import env
import pdb

##def WalkErr(onerror):
##    try:
##        import errno
##        print errno
##        if errno!=errno.ENOENT: raise
##
##        elif errno == errno.EACCES:
##            print filename,"cannot open"
##        else:
##            print errno
##            print filename, "Not Found....Continuing"
##    except:
##        print "Something went wrong"

WorkspaceRoot="W:/GIS_Scratch"

env.workspace=WorkspaceRoot

#List workspaces in root directory
workspaces=arcpy.ListWorkspaces ("*","Folder")
for workspace in workspaces:
    env.workspace=workspace
    print workspace

#Find workspaces named Drawings
    workspaces1=arcpy.ListWorkspaces("*rawing*","Folder")
    print workspaces1

    for workspaceNext in workspaces1:
        arcpy.env.workspace = workspaceNext
        print workspaceNext

#start stepping through code
#    pdb.set_trace()

# Walk through directory to find CAD polylines (Note: "CadDrawing" datatype does not return any polylines, so had to use "FeatureClass")
        try:
            for dirpath, dirnames, filenames in arcpy.da.Walk(workspaceNext,topdown=True, datatype='FeatureClass',onerror=None,type='Polyline'):
                try:
                    print dirpath
                    for filename in filenames:

                        try:

                            arcpy.env.workspace=filename
                            desc = arcpy.Describe(dirpath)
                            featureTest=(os.path.join(dirpath,filename))
                            print filename
                            print featureTest
                            print desc.extent.XMin
                            filename.close

                        except arcpy.ExecuteError:
                            arcpy.AddError(arcpy.GetMessages())
                        except:
                            arcpy.AddError("Non-tool error occurred")
                except OSError as e:
                    print e.errno
                    print e.filename
                    print e.strerror
                    continue

                except Exception as e:
                    print "An error has occurred"
                    print e
                    continue

        except OSError as e:
            print e.errno
            print e.filename
            print e.strerror
            continue
        except Exception as e:
            print "An error has occurred"
            print e
            continue


texts: QGIS Snapping is not aligned with layers

texts:
I was trying to draw lines based on hex grid, but snapping give me wrong position of hex corner.

As you can see, I used snapping-advanced-25 pixel from hex grid, but it is not aligned hex grid. It is off from hex grids a bit.

How can I solve this issue?


texts: Reduce image collection to get daily sum from hourly precipitation GSMaP 12:00h pm to 11:00h am

texts:
I am trying to adapt a script made available by colleagues from the GEE of the GSMaP image collection. The script obtains data (as a table) of the points of each image available in a given period and adds the hourly data per day. In my case, I want to compare this data (daily sum) with rain data that is collected from 12:00pm to 11:00am (Brazil). Therefore, you would need help for the script to add, for example, from 12:00pm 01/01/2019 to 11:00am 01/02/2019. Please find the script that needs to be improved below:

var p1 = /* color: #d63000 */ee.Geometry.Point([69.15, 54.8331]),
    p2 = /* color: #98ff00 */ee.Geometry.Point([70.917, 54.433]),
    p3 = /* color: #0b4a8b */ee.Geometry.Point([66.967, 54.367]);

// Collection of point
var pts = ee.FeatureCollection(ee.List([ee.Feature(p1),ee.Feature(p2), ee.Feature(p3)]))

// Start and End Dates
var inidate = ee.Date.fromYMD(2014,1,1)
var enddate = ee.Date.fromYMD(2014,3,1)

// Difference between start and end in days
var difdate = enddate.difference(inidate, 'day')

// Import GSMaP data
var gsmap = ee.ImageCollection('JAXA/GPM_L3/GSMaP/v6/reanalysis')
.filterDate(inidate, enddate)
.select('hourlyPrecipRateGC');

// Time lapse
var lapse = ee.List.sequence(0, difdate.subtract(1))
var inidate = ee.Date('2014-01-01')
var listdates = lapse.map(function(day){
  return inidate.advance(day, 'day')
})

// Iterate over the list of dates
var newft = ee.FeatureCollection(listdates.iterate(function(day, ft) {
  // Cast
  ft = ee.FeatureCollection(ft)
  day = ee.Date(day)

  // Filter the collection in one day (24 images)
  var day_collection = gsmap.filterDate(day, day.advance(1, 'day'))

  // Get the sum of all 24 images into one Image
  var sum = ee.Image(day_collection.sum())

  // map over the ft to set a property over each feature
  var temp_ft = ft.map(function(feat){
    var geom = feat.geometry()
    var value = sum.reduceRegion(ee.Reducer.first(), geom, 30).get('hourlyPrecipRateGC')
    return feat.set(day.format('YYYY-MM-dd'), value)
  })

  // Return the FeatureCollection with the new properties set
  return temp_ft
}, pts))

print(newft)

// Export
Export.table.toDrive(newft,
"KZ_GSMaP_14-15",
"KZ_GSMaP_point data",
"KZ_GSMaP_14-15")


texts:
I am using r.watershed to get a stream map (which I assume is rivers). However the output rivers is always 1 pixel wide. I would've assumed that if I run this on say Iraq I would see the Tigris and Euphrates stand out / be abit larger.
Is there a way to tweak r.watershed to do this or is there another GRASS tool for this?
I am running on SRTM elevation maps if that makes a difference.


texts: Error saving geodaframe to shp: AttributeError: 'Series' object has no attribute 'to_file'

texts:
I have been struggling with a problem related to Getting AttributeError: 'DataFrame' object has no attribute 'to_file' from GeoPandas even though file converted into GeoDataframe instead of Dataframe but it is not exactly the same.
I created a point geopandas dataframe and saved it to a shapefile without a problem
pnt_data = {'geometry': [<shapely.geometry.point.Point object at 0x7fa9820ffd10>, <shapely.geometry.point.Point object at 0x7fa95739ea90>,...

gp_df = gpd.GeoDataFrame(pnt_data, crs=crs)
gp_df.to_file(path_to_file, driver="ESRI Shapefile")

Now I want to save a line shapefile with the same points. Based on what I've read, I tried the following:
lineCoords = [(500191.6, 7656332.248), (500222.84883889696,
               7656358.309988857),...]
df = gpd.GeoDataFrame(lineCoords, crs=crs)

which got me the following dataframe
                                 geometry
0                  (500191.6, 7656332.248)
1  (500222.84883889696, 7656358.309988857)
2   (500238.1988799781, 7656371.130918575)
3  (500254.99698246736, 7656381.985583601)
4    (500273.2863161564, 7656390.07880751)
5  (500292.37426232523, 7656396.049595658)
6  (500312.1612412515, 7656398.9608608475)
7  (500331.89439794194, 7656402.217013037)
8  (500351.87729180406, 7656403.044027512)
9  (500371.78081168584, 7656401.081912234)

The only attribute that interests me are the coordinates and they are there. However, when I tried to save this to a shapefile using the same procedure I used for points
gp_df.to_file(path_to_file, driver="ESRI Shapefile")

and got the error message
AttributeError: 'Series' object has no attribute 'has_z'

Based on a similar problem I found that the solution could have been
gp_df['geometry'].to_file(filename, driver="ESRI Shapefile")

At this point I got the error in the title:
AttributeError: 'Series' object has no attribute 'to_file'

I also tried to convert the list of tuples with the coordinates to a linestring and do the same procedure
finalData = LineString (lineCoords)
gp_df = gpd.GeoDataFrame(finalData, crs=crs)

In this case I got
ValueError: DataFrame constructor not properly called!

Some recommended to use the raw list instead of the linestring, which was why I started the path above.


texts: Adding time stamps to animation slides with OSS tools

texts:
I would like to create an animation from a set of PNG files and apply labels (time lapse values) on each of the PNG files. The PNG files were actually heat maps that were generated by interpolation.
I have created the heat maps in 2 steps:

I Created a QGIS a vector layer that stores timestamp and temperature value pairs for several points on the map. I also added in the vector layer a time stamp for each data point. The header looks something like:

point1, lat, long, time1, temp11, time2, temp12, ......
point2, lat, long, time1, temp12, time2, temp22, ......

I interpolated each vector layer point for each time stamp (time 1, time2 etc) in GRASS using v.surf.icw that resulted in creating one raster layer for each time stamp, that I exported as a PNG files.

I combined all the created images in a movie with the GRASS 7 GUI animator. Everything looks good except that I would like to have the time lapse label printed on each slide/frame.
In order to add the labels on the slides, I would like now to create a script that reads the time values from the vector layer that was detailed above and write/print the selected values then on the PNG images.
One way I found is to use d.text in grass as documented here. However, this example only shows how to create new slides but  not how to write on existing ones. In addition, the time stamp information needs to be extracted/read from the vector layer. I assume with d.vect or something similar may work.
Therefore my question about how put a label on a animation slide has two aspects:

Extract an attribute value from a vector layer?
How to write the value of the attribute with d.text on an exiting PNG image?


texts: Extract exact boundary of DEM from another raster

texts: I have a DEM and I want to extract this DEM by using another raster. Basically I need only those DEM pixels that are covered by second raster. When I try to clip, the output gives the whole extent of the second raster. Extract by mask is also not working.


texts: Identification of consecutive points within a specified buffer

texts:
I have a point file, hourly relocations of an animal, and I want to be able to place a buffer around each point and calculate the number of subsequent points which are within the buffer zone. I am looking for a method which will work along the point file, like a moving window, which will only count those subsequent points which are within the buffer zone.

For instance at point 10 I place a buffer of 1500m and I want to know whether point 11 is within the buffer and if so then whether point 12 is within the buffer and so forth. I don't want to know whether point 52 is within the buffer zone unless all of the previous points have been within the buffer. I also don't want to know whether points 9 or 8 etc. are within the buffer.

I have found and tried an additional toolbox called "moving window point analysis" which works as a moving window on the point file. This works well, but very slowly, and includes all points that are within the buffer zone even if they are not consecutive points. I can't find a way to just make it look at consecutive points.

I would like a method which will provide an output table as I have a lot of data points to look at in this way.

I am using ArcGIS 10. Any help that anyone can provide would be greatly appreciated.


texts: Getting python geoservice ExportToJPEG to return image as job result?

texts:
I have an asynchronous geoprocessing service written in python that takes a layer from a mapservice and generates a jpeg image out of it.

To achieve this I am using ExportToJPEG

I need said geoprocessing service to return the image as a job result. How can I achieve that?

import arcpy
q ={
"mapOptions": {
"extent": {
            "xmin" : -79.50,
            "ymin" : -4.60,
            "xmax" : -66.50,
            "ymax" : 12.8,
            "spatialReference" : {
                "wkid" : 4686
            }
        }
},
"operationalLayers": [{ "url":"....."}]
}
result = arcpy.mapping.ConvertWebMapToMapDocument(q);
arcpy.mapping.ExportToJPEG(result.mapDocument, path_to_image)
del result


texts: How or where do I find the date of an NAIP orthophoto, specifically a California NAIP 2016? I need to know when an image was created.


texts: Converting {"type":"MultiLineString","coordinates":[[[x1,y1],[x2,y2]]]} to geography column using PostGIS?

texts:
I have a Postgresql table with a text column that has data like the below.

{"type":"MultiLineString","coordinates":[[[-122.995882414,44.05766941700006],[-122.99818210099994,44.05768878400005]]]}

I have added a new column to the table with type Geography and I would like to move this multilinestring data to the new geography column.

I need to find the centroid of the mulitilinestring in the form of lat/lng so that distance can be calculated from a give address. The table has around 139k rows.


texts: Why does QGIS crash when my PostGIS table has a primary key?

texts:
Configuration:

PostGIS Library: 2.0.3 GEOS: 3.3.8-CAPI-1.7.8 Proj: Rel. 4.8.0, 6 March 2012

PostgreSQL 9.2.2, compiled by VC++ build 1600, 64 bit, running on Win2012 server 64 bit

QGIS 2.2.0 - 32 bit Windows, running on Windows 7 64 bit.

---> The tl;dr

Another way of looking at this issue is that any primary key added to the spatial table causes QGIS to not find a primary key on the table.

<---

I've been struggling with a QGIS/PostGIS related problem for a few days now.

My PostGIS database has a spatial table with integer column called "gid" that is assigned a primary key constraint. That spatial table can be loaded as a spatial layer in QGIS, but opening its table causes QGIS to hang and the table doesn't display. Using the select tool to select one feature within the map display will select all features.

I've also noticed in the layer properties that QGIS is unable to identify a primary key.

I found that if I remove the primary key constraint from the column (now the table does not have any primary key) the spatial table opens fine in QGIS and the issue is no longer apparent, but other issues occur because QGIS cannot identify a primary key for the table.

I've found a work around - if I create a numeric column called oid and populate it with a sequence and do not assign a primary key to the table, QGIS can identify the oid column as a primary key. It doesn't seem to identify a similarly built column name of "gid" as a primary key, however.

I've tried a few different configurations of PostGIS and QGIS.

I've used QGIS 2.0 - 64 bit for Windows and experienced the same issue. I've tried another database server with identical configurations, functions, and extensions, and connecting QGIS to that server and loading spatial layers does not exhibit the issue; QGIS is able to identify a primary key column.

Has anyone experienced issues with QGIS being unable to identify a primary key column when the table has a primary key in Postgres? I've searched through stack overflow and found a few issues that seem similar, but not identical to this problem.

This issue appears similar, but I've used the DB manager and not experienced any success in this situation. When loading a spatial layer from PostGIS via "Add PostGIS Layer" QGIS doesn't provide a drop down box to select a primary key for any spatial layers.

https://hub.qgis.org/issues/9803


texts: Determine PostGIS geography polygon orientation

texts:
How can you determine the orientation of a PostGIS geography polygon?  I'm looking/hoping for a function that would tell me which side is the interior (left or right) as you traverse the vertices of the outer shell in order.

I have seen the ST_Orientation function, but it only works for geometry.  I have also seen the ST_ForceRHR and ST_ForceLHR functions.  But again, they only work for geometry.

I have data that wraps around the poles, so naive casting to geometry won't work for me.  (Or at least I don't think so.  Maybe this isn't true.)

I'm a reasonably good programmer, so if no such function exists then an algorithm outline would be fine too.


texts:
I followed the instructions how to build a dam by whuber (How to build an artificial 'dam' across a river?)

"draw a line segment"

"compute the Euclidean distance to this line segment"

"divide that by 5" for the correct slope (no picture here, looks like the result of step before divided by 5...)
"subtract that from the elevation of the top of the dam" (I want my dam to be 2 m high)

"combine it with the original DEM"

My result (damdtm)

It is correct for the top of the dam:

But incorrect for the rest of the processing extent:

Shouldn't there be a step before step 5 where you cut the dams footprint from the raster before combining it with the original DEM?


texts: Increasing extremely small text in ArcMap?

texts:
I'm taking an introductory GIS course and the university provided us with one year student licenses so that we can work from home. The issue is the text in the table of contents and the toolbar(?) (File, Edit, View, Etc.) up top are in very small font; so are the icons and text in the table of contents.

For a while a lot of the toolbar icons were also incredibly small but I figured that out by using the large icon sizes, however the icons are still small within a tool.

A few months ago I had the same issue with QGIS. I originally deleted the program but since then have redownloaded QGIS and no longer have that issue. To me, this hints at the issue being a laptop problem? I don't know though, just a guess. I'm running windows 10 on an Asus zenbook.

UPDATE:

Adjusting screen resolution has helped but some icons and text (like in arctoolbox) are still really small:

I tried referencing these two posts but am still  having trouble:

Size of icons and texts in ArcMap?
Where does ArcMap GUI get font size and docking guide positions from?


texts:
Right I've got this JS function

   function AreaCompare(sql, UIDfield, Geogname, Table, sourceName, sourceID, destName, destID){

    CompareObj = {};
    CompareObj.sName = sourceName;
    CompareObj.sID = sourceID;
    CompareObj.dName = destName;
    CompareObj.dID =destID;

    sql.execute("select ST_Area(the_geom::geography)/1000000 as SA FROM "+Table+" where "+UIDfield+"="+sourceID)
    .done(function(data) {

        CompareObj.sArea=JSON.stringify(data.rows[0].sa);

        sql.execute("select ST_Area(the_geom::geography)/1000000 as DA FROM "+Table+" where "+UIDfield+"="+destID)
        .done(function(data) {

            CompareObj.dArea=JSON.stringify(data.rows[0].da);
            alert(JSON.stringify(CompareObj));
            return CompareObj;

        })
    })
}

It appears to work correctly and the alert box flashes up with a readable version of the object that is as expected.

however when I examine the returned object it is undefined, essentially empty/null.

What trick am I missing a trick here? It originally thought it was something related to scope and the done function.

I did allot of research and I realised I needed to nest the CartoDB SQL statements, but that didn't fix the problem.

Any suggestions?


texts: I'm using Debian 7 and yesterday i had an update of Qgis 2.8. Since then I can't change font size and type, at labeling and print composer, whatever i pick it returns to the default  "Free sans" size 11. Any suggestions?


texts: Cannot create pgr_PointsAsPolygon from pgr_DrivingDistance function?

texts:
I have successfully build a pgr_DrivingDistance query,

select
    id::int4,
    lon::float8 as x,
    lat::float8 as y
from
pgr_DrivingDistance(
    'select gid as id, source, target, length as cost from ways',
    23199,
    0.01,
    false
) as isokron inner join ways_vertices_pgr on (isokron.node = ways_vertices_pgr.id)

to have successfully have an id, x, and y column for the output.

I want to use that output into pgr_PointsAsPolygon function, which will create a polygon from the (x, y) pair above.

But I'm having difficulty in building the combined. I try these, but without result.

select pgr_pointsAsPolygon(
'select
    id,
    lon as x,
    lat as y
from
    pgr_drivingDistance(
        ''select gid as id, source, target, length as cost from ways'',
        23199,
        0.01,
        false
    ) as isokron inner join ways_vertices_pgr on (isokron.node = ways_vertices_pgr.id)');

It keeps telling me this error:

ERROR:  syntax error at or near "select"
LINE 7:    'select gid as id, source, target, length as cost from wa...
            ^
QUERY:  SELECT x, y FROM pgr_alphashape('select
    id,
    lon as x,
    lat as y
from
    pgr_drivingDistance(
        'select gid as id, source, target, length as cost from ways',
        23199,
        0.01,
        false
    ) as isokron inner join ways_vertices_pgr on (isokron.node = ways_vertices_pgr.id)', 0)
CONTEXT:  PL/pgSQL function pgr_pointsaspolygon(character varying,double precision) line 17 at FOR over EXECUTE statement
********** Error **********

ERROR: syntax error at or near "select"
SQL state: 42601
Context: PL/pgSQL function pgr_pointsaspolygon(character varying,double precision) line 17 at FOR over EXECUTE statement

I've tried create the polygons (pgr_pointsAsPolygon) using all the vertices (without using pgr_drivingDistance) and it worked perfectly!

Whats wrong with the code?
Is the double single-quote escape wrong?


texts: Overpass API python wrapper to run requests on the endpoint of overpass-turbo

texts:
we have with Overpass API python wrapper a thin Python wrapper around the OpenStreetMap Overpass API  https://github.com/mvexel/overpass-api-python-wrapper

we have some Simple example:

import overpass
api = overpass.API()
response = api.Get('node["name"="Salt Lake City"]')

Note that you don't have to include any of the output meta statements. The wrapper will, well, wrap those.

we will get our result as a dictionary, which represents the JSON output you would get from the Overpass API directly.

print [(feature['tags']['name'], feature['id']) for feature in response['elements']]
[(u'Salt Lake City', 150935219), (u'Salt Lake City', 585370637), (u'Salt Lake City', 1615721573)]

we an specify the format of the response. By default, we will get GeoJSON using the responseformat parameter. Alternatives are plain JSON (json) and OSM XML (xml), as ouput directly by the Overpass API.

response = api.Get('node["name"="Salt Lake City"]', responseformat="xml")

updated question: can we also get cvs - can we perform a request like below with the python wrapper to the endpoint of overpass turbo?

[out:csv(::id,::type,"name","addr:postcode","addr:city",
"addr:street","addr:housenumber","website"," contact:email=*")][timeout:30];
area[name="Madrid"]->.a;
( node(area.a)[amenity=hospital];
  way(area.a)[amenity=hospital];
  rel(area.a)[amenity=hospital];);
out;


texts: QGIS 2.18.2 on mac cannot style or label points from csv vector layer

texts: Suddenly with QGIS 2.18.2 on Mac OSX 10.12.3, when I import a CSV as a vector layer, the layer properties > style and >labels dialogs are completely blank. It was working before. Same files work fine with QGIS 2.14.0 on windows 10. what gives? Is there a setting I can tweak? a dependency? like JavaScript or a plugin?


texts:
I'm working with ArcGIS for desktop. I have received data from a third party and am trying to update my CAD system via GISLink. In order to connect the new data with the currently-used data, certain fields are required, the most important of which being a unique External Key field. In the new data tables, this field contains NULL values or is left out completely. In fact, there are no unique identifiers that match up between tables. I have so far been unable to reach the person who provided the data.
I am new to using GISLink and data management in general. Is there a way of working around this short of manually matching and calculating the leftovers? There are tens of thousands of rows, so this is not really a viable option.
It seems to me there is nothing I can do until I can contact the source, but I wanted to make sure I'm not missing something obvious.


texts: Automatically refresh my wfs layer in Leaflet when adding data in Postgis

texts:
I’ve built a web-sig based on postgis, geoserver and leaflet.
My polylines, …, are stored in Postgis and I call them through Geoserver in Leaflet as wfs layers like :

var defaultParameters = {
    service : 'WFS',
    version : '1.0.0',
    request : 'GetFeature',
    typeName : 'database:table’,
    outputFormat : 'text/javascript',
    format_options : 'callback:callEP',
    srsName : 'EPSG:4326',
    };

var parameters = L.Util.extend(defaultParameters);

var ajaxEP = $.ajax({
    url : owsrootUrlAssainissement + L.Util.getParamString(parameters),
    dataType : 'jsonp',
    jsonpCallback: 'callEP',
    success : function (layerEP) {
        EP = L.geoJson(layerEP, {
            style: function (feature) {
                return {
                    color: "lime",
                    fillColor: 'lime',
                    opacity: 1,
                    weight: 3,
                };
            },
     });
   }
});

I’m also using leaflet draw to draw new polylines and I send them into Postgis using ajax and php :

map.on('draw:created', function (e) {
        layerEP = e.layer;
        var coordinatesEP = layerEP.getLatLngs();
        var sql_coordinatesEP =[];
        for (var i in coordinatesEP){
                sql_coordinatesEP.push(coordinatesEP[i].lng+" "+coordinatesEP[i].lat);
            }

        dialog_create_EP.dialog("open");

        $("#formulaire").submit(function(){
                $('#coordsCreateEP').val(sql_coordinatesEP);
                $('#longueurCreateEP').val(longueurEP);

                                    var formData = {
                                      'type' : $('select[name=type]').val(),
                                      'diametre' : $('select[name=diametre]').val(),
                                      'coordsCreateEP' : $('#coordsCreateEP').val(),
                                    };
                                    $.ajax({
                                      url:'create_EP.php',
                                      type: 'POST',
                                      data: formData,
                                      success : function(data, response) {
                                        dialog_create_EP.dialog("close");
                                    }
                               });
                         });
    });
});

php :

$type = $_POST['type'];
$diametre = $_POST['diametre'];
$coordsCreateEP = $_POST['coordsCreateEP'];

try
{
   $bdd = new PDO("pgsql:host=$host;dbname=$database",$user,$pass);
   $bdd->setAttribute(PDO::ATTR_ERRMODE, PDO::ERRMODE_EXCEPTION);

$req = $bdd->prepare('INSERT INTO database.table(type,diametre,coords) VALUES(?,?,?)');

$req->execute(array($type,$diametre,$coordsCreateEP));

   }
catch(Exception $e)
{
    echo "Connection a la BDD impossible : ", $e->getMessage();
 die();
}

Everything works fine, and when I’m sending new polylines via ajax, my page is reloaded and my new polylines appears when I’m adding their layer.

Ok, now, What I would like to do would be to not refresh all my page when sending new polylines but just the layer with the new polylines. So is there a way to refresh a leaflet layer automatically when adding data into Postgis without reloading all my page ?


texts:
I am trying to merge polygons using shapely. Here is the piece of code that is being used. This code, fully functional earlier, doesn't seem to work (probably because my python packages were updated).

#open shapefile
fshp = fiona.open(shpFileWetlandPath)
shapesTotal = []
#create collection of shapes to merge
for counterShp in range(0, len(fshp)):
    shapefile_record = fshp[counterShp]
    shape = shapely.geometry.asShape(shapefile_record['geometry'])
    print("adding shape " + str(counterShp))
    shapesTotal.append(shape)
print("shapes added. attempt union")
shapeMerged = cascaded_union(shapesTotal)

This fails with the error message:

Assertion failed: (!"should never be reached"), function itemsTree, file AbstractSTRtree.cpp, line 372.

I am using shapely 1.5.16 and python 2.7.11.


texts: How to ensure transparent background when exporting a map from ArcGIS Pro to PDF?

texts:
Is there a way to export a PDF from ArcGIS Pro 2.7.2 to PDF, ensuring there truly is no background? When I export from ArcGIS Pro to PDF, even when I specify that the map frame be transparent and the layout be transparent, when I export the map to PDF it still maintains a white background layer.
To see this in action, I went to Adobe Illustrator. I opened up 1 map PDF where I did include a white background, and then opened another map PDF where I specified a transparent/no background. I made a teal rectangle in Illustrator and dragged it to the bottom of each map to see what was covering it, pictures below.
I couldn't find much online, only this link: https://community.esri.com/t5/cartography-and-maps-questions/how-to-make-a-layout-background-transparent-in-arcgis-pro/td-p/774159
1st pic, in Pro I specific this map have a white background. It has 2 white backgrounds layers in AI, highlighted.

2nd pic, in Pro I specified the map frame have no background. It has 1 white background layer in AI, highlighted.


texts: how to change the mouse cursor into pan when i select the pan icon in openLayers?

texts:
I want to change the mouse cursor into a hand icon when I select to pan.
I have tried this code:

map.addControl(new OpenLayers.Control.PanZoomBar());
map.addControl(new OpenLayers.Control.MouseToolbar());
map.addControl(new OpenLayers.Control.LayerSwitcher({'ascending':false}));
map.addControl(new OpenLayers.Control.Permalink());
map.addControl(new OpenLayers.Control.Permalink('permalink'));
map.addControl(new OpenLayers.Control.MousePosition());
map.addControl(new OpenLayers.Control.OverviewMap());
map.addControl(new OpenLayers.Control.KeyboardDefaults());
map.addControl(new OpenLayers.Control.Dragpan());
map.addControl(new OpenLayers.Control.Pan());
map.addControl(new OpenLayers.Control.PanPanel());

OpenLayers.Control.DragPan = OpenLayers.Class(OpenLayers.Control, {
type: OpenLayers.Control.TYPE_TOOL,
panned: false,
interval: 25,
documentDrag: false,
draw: function() {

this.handler = new OpenLayers.Handler.Drag(this, {
                "move": this.panMap,
                "done": this.panMapDone
            }, {
                interval: this.interval,
                documentDrag: this.documentDrag
            }
        );
    },
panMap: function(xy) {
        this.panned = true;
        this.map.pan(
            this.handler.last.x - xy.x,
            this.handler.last.y - xy.y,
            {dragging: this.handler.dragging, animate: false}
        );
    },
 panMapDone: function(xy) {
        if(this.panned) {
            this.panMap(xy);
            this.panned = false;
        }
    },

    CLASS_NAME: "OpenLayers.Control.DragPan"
});
map.addControl(control);


texts: Creating new raster based on spatial extent and classes of one raster and average of another using ArcMap?

texts:
I have a raster that is split into classes and I wanted to make a new raster with this raster's spatial extent and classes but displaying the average values from a different raster.

How do I do this?


texts: MBG_Orientation: ranges from 0-180, instead of 0-360? (ArcGIS, Minimal Bounding Geometry)

texts:
I'm trying to understand the concept of the Orientation of the Minimal Bounding Geometry: Minimum Bounding Geometry (Data Management toolbox, ArcGIS Desktop 10.1).

I run the tool Minimal Bounding geometry -> Rectangle_by_width -> "Add geometry characteristics as attributes to output". For every polygon I've got back: MBG_Width, MBG_Length, MBG_Orientation. As it is written in the tool description: "the orientation angles are in decimal degrees clockwise from north". My resulting angles range, however, from 0-180, not to 360, as I would expect.

How can I interpret my results? I doubt to have only eastern orientations in my data. Is there any simple drawing explaining how the MBG_Orientation is calculated in ArcGIS, i.e. from which point, to which direction? Or why the resulting orientation is from 0-180?


texts:
I'm building dynamic cartographic symbols to display inside ArcMap, but they are displaying "wrong".

Is there a guide somewhere that will show what raster op code does? I mean, not only show, but perhaps demonstrate the symbology. The esri pages has confusing descriptions like "it's the inverse of rasterOpCodeXOR" for example.

I've always struggled with this. Any change there are good resources on this?


texts: area weighted average for overlapping polygons

texts:
Possible Duplicate:
  How to use NASS county level crop data to estimate crop acreage in hydrologic regions?

It seems like an easy problem but I don't know how to go about it.
I have two layers for British Columbia, one divides the province into 63 polygons and the other layer divides the province into 20 polygons. What I want to do is assign the values from the the 20 polygons to the 63 polygons. Since the boundaries are not the same I would like to assign an area weighted average, so that the 63 polygons have the values from the 20 polygons depending on how much of their area overlaps. Does something like that exist or what other approach is usually applied in that case? I feel like trying to combine data with different spatial extent is a common task?


texts:
I have tiles generated by Maperitive from an OSM file, and there's a nice tiles.json generated with them which specifies neatly the tiles-path, bounding box, minzoom, maxzoom, center, etc. The format is identical to the tileJSON specification given here: https://www.mapbox.com/guides/an-open-platform/#tilejson

But it just blows that the example stopped a step short and inline'd the JSON code into the example's source code, hence rendering it (pun intended) useless. No clue on https://github.com/mapbox/tilejson-spec either. All web searches for this lead to these sites; the rest digress into GeoJSON which is something else entirely. Why in the world would anyone go through the trouble of using JSON if one wasn't importing it from an external file? (no, please don't try to answer that.. it's not the main question. Focus!)

So, a humble request: Please re-write this example code given on mapbox.com, to show an external JSON file being loaded. (And please put in a word with someone there if you know them to post a proper example on the website)

<!DOCTYPE html>
<html>
<head>
  <meta charset=utf-8 />
  <title></title>
  <script src='https://api.tiles.mapbox.com/mapbox.js/v2.1.4/mapbox.js'></script>
  <link href='https://api.tiles.mapbox.com/mapbox.js/v2.1.4/mapbox.css' rel='stylesheet' />
  <style>
    body { margin:0; padding:0; }
    .map { position:absolute; top:0; bottom:0; width:100%; }
  </style>
</head>
<body>
<div id='map' class='map'> </div>
<script>
L.mapbox.accessToken = '<your access token here>';
var tilejson = {
  "tiles": [ "https://api.tiles.mapbox.com/v4/examples.map-i86l3621/{z}/{x}/{y}.png?access_token=pk.eyJ1IjoidHJpc3RlbiIsImEiOiJuZ2E5MG5BIn0.39lpfFC5Nxyqck1qbTNquQ" ],
  "minzoom": 0,
  "maxzoom": 18
}
L.mapbox.map('map', tilejson, {
    scrollWheelZoom: false
}).setView([43.6519,-79.3852], 15);
</script>
</body>
</html>

Edit: Included jquery.min.js and changed code as advised; am getting a "not well-formed" error in console. I reckon jquery doesn't like tileJSON. This is why full examples ought to be given :P. Please see and advise what needs to change here. Is there a more generic JSON-loader we can use instead of jquery? I believe we're very close to getting this.

var tilejson;
$.getJSON('pune/tiles.json', function(data) {
  tilejson = data;
});
console.log(tilejson.name);

L.mapbox.map('map', tilejson, {
    scrollWheelZoom: false
});

Here is the content of tiles.json (loading local tiles, generated by Maperitive; tested to be valid JSON as per http://codebeautify.org/jsonviewer):

{
    "tilejson": "1.0.0",
    "name": "Pune 1 to 16",
    "description": "Made with Maperitive",
    "attribution": "Map data © OpenStreetMap contributors",
    "tiles": [
        "pune/{z}/{x}/{y}.png"
    ],
    "minzoom": 1,
    "maxzoom": 16,
    "bounds": [
        73.59366818720622,
        18.33281593564951,
        74.07630689997964,
        18.77938751279949
    ],
    "center": [
        73.83498754359293,
        18.5561017242245,
        1
    ]
}


texts: Merge with ogr2ogr creating multiple feature classes despite consistent -nln flag

texts:
I have a Python script that creates hundreds of individual shapefiles and in the end, merges them iteratively to one big geodatabase using ogr2ogr.

I want one geodatabase output with a single feature class and all of the polygons from the input shapefiles in that feature class, but that is not what is happening and I cannot figure out why.

This is what I am doing:

# if output.gdb does not exist:
ogr2ogr -nln output -a_srs EPSG:4326 -t_srs EPSG:4326 -f "FileGDB" output.gdb input1.shp

# if output.gdb does exist (i.e. not on the first iteration):
ogr2ogr -nln output -a_srs EPSG:4326 -t_srs EPSG:4326 -update -append -f "FileGDB" output.gdb input2.shp

# And repeat the second command for each inputX.shp (in my case, 109 shapefiles, 66959 total points)

Please note that the names have been changed for readability. The input shapefiles are not numbered like that, so the feature class names in the resulting output.gdb are not a result of the input names.

This is what the output (with the names again changed to be more readable and consistent with the above) looks like when I run ogrinfo on output.gdb [broken up into 2 screenshots]:

When viewing the individual feature classes in Arc, it appears that each represents one of the inputs. What am I doing wrong here? What is with all the feature classes, their names, and why are there only 100 of them as opposed to 109? How can I rectify this to get my desired end result?

UPDATE:

I wanted to check and be sure the field names were consistent throughout. I compared the field names from 2 random input shapefiles, and they were the same (as I would expect for all of them). Then I checked against the output geodatabase and discovered all of my columns from the shapefiles that began with a number were different in the gdb (they all had been changed to have an underscore before them). I found out the gdb's cannot have field names that begin with a number. Could this be why the above commands did not want to merge the shapefiles into the same feature class within the output geodatabase? I will run a test to see.

UPDATE 2:

I ran 2 tests: one merging the 3 shapefiles with some columns beginning with a number, and another merging the 3 shapefiles where those columns had been changed to begin with a letter. The resulting gdb of the latter test had field names that matched the input shapefiles exactly, but both tests resulted in the same issue where the output gdb had 3 feature classes instead of the results being in one single feature class. So that was not the problem.

UPDATE 3:

Tried a bunch of suggestions from this answer, and still same issue. I also tried ogrmerge.py and I can't write to a gdb in update mode. At a loss.


texts: pg_tileserv database connection string specify port

texts:
I am trying to connect a DB to pg_tileserv following these directions https://github.com/CrunchyData/pg_tileserv/blob/master/README.md
it says after download for windows:
SET DATABASE_URL=postgresql://username:password@host/dbname
pg_tileserv.exe

but there seems to be no option to pass in the port. I have to connect to another port that is not 5432


texts: How could I create a heatmap with carto.js 4.2.X, is it possible?

texts: I am trying to create a heatmap with the latest version of carto.js and gmaps, but I have not succeeded. I don't know if it is possible with this new version. With cartoDB.js it was simpler, because the layer type was indicated in TORQUE.


texts: Analyzing spatial pattern in time using ArcMap?

texts:
I am analyzing a dispersion of unemployement rate at municipality level in Czechia between 2002 and 2016.
I have chosen LISA analysis which identifies high-high and low-low clusters.

Now I would like to identify and measure whether the municipalities (polygons) are in the same clusters for the whole period. Which municipalites were in high-high cluster at the beginning and have changed the position to low-low clusters after years? I would like to identify which municipalites are stable and are stucked in the one cluster and which ones have moved to any other clusters.

Unfortunately, I do not know how to do it. My only idea is to join all LISA results into one shapefile, then export it to the Excel and then try to analyze it there using excel functions like "if". I believe, there must be a better (and statisticaly correct) method in ArcMap...

Sorry for my english and my poor knowledge of GIS.


texts: Reconciling elevation difference between SRTM and old US Army map in Laos?

texts:
In QGIS I georeferenced a 1970 US Army contour map for part of southern Laos (20 m intervals) with the following indicated properties:

Horizontal datum Indian datum 1960
Spheroid=Everest
Verical datum ="approximate mean sea level"

It shows elevations about 12 m lower than my SRTM tif which I set to WGS84 EPSG:4326. The Army map fits other information I have for my project better. I guess I could just adjust the base for my project by 12 m but I hate to do this without understanding the big difference.

I'm new to all this. Any ideas on how to reconcile the two?


texts: Doing automatic vectorization in ArcMap?

texts:
I'm using ArcMap 10.6; ArcScan. I am trying to run an automatic vectorization on my .tif file. I've cleaned up my raster, and am doing Start Editing - Vectorization - Show Preview (looks good) - Generate Features ... and this is where it has the progress bar, like it's creating them, and my map 'refreshes', yet... nothing is in my line layer.

I have both a parcelline and parcelpolygon layer, both are file geodatabase feature classes. They are what I select when Generating Features. When I Start Editing, I'm selecting my file geodatabase that contains those two layers.

I've tried all the trouble shooting I know -- I even tried creating a shp file for the lines, like another post recommended, but to no avail. Help?

I believe my tif falls in the same spatial domain. This was the previous post Error "generate features" with ArcScan vectorization mentioned.


texts: I would like to take excel tables and place them in a label or annotation so that I no longer have to copy and paste them into the figure layout as picture objects. I know this isn't possible in ArcMap but is it possible in QGIS?


texts: Tiling failing near the end of generating base level tiles


texts: How to find all intersections and count lines arriving at each intersection?

texts:
I'm having a hard time getting off the ground with this, and am trying to move away from ArcGIS / arcpy.

I have a shapefile with a network of polylines.  I want to find/count all 2-way intersections, 3-way intersections, 4-way, etc.  This would involve, I think:

splitting all lines at intersections
deleting identical segments (with some spatial tolerance)
create points at the end of each new line segment
buffer the line ends, and count the number of points near/at each line end

I would like to do this in a python script. Is this overall plan easily implementable in pyQGIS? An outline of the calls I would be using would be very helpful.


texts: How to create histogram to visualize feature's value in QGIS 3?

texts: Is there any tool or plugin can create a plot like this (generated by Arcmap 10.7)? My QGIS software's version is 3.10.  What I know is that the QGIS built-in tool, basic statisc for field, provides only some statistics, no plot available.


texts: How to create multipe labels on polygons?

texts:
Is there any tool or function to create muliple labels on polygons?

Edit:

Meanwhile I've created a point layer for label placement using the refFunctions plugin.

geomwithin('polygon', 'label')

Unfortunately QGIS gets quite slow having a few hundred point features.


texts: QGIS Invalid Geometry when importing a CSV file

texts:
I've been learning how to use QGIS for the past couple of weeks, but am having some trouble importing some of my own data as a CSV file. My data looks like:

Which is giving me the following error:

  49 records discarded due to invalid geometry definitions.  The
  following lines were not loaded into QGIS
  Invalid X or Y fields at line 2 (and so on)

I'm using QGIS on Linux but originally saved the data in Windows/Excel and thought that may have something to do with it, so tried saving it in Libre but no luck there. I also read on a similar question on here that foreign characters in other fields was the problem, so I removed (äüß) them but no luck.

Edit 2:


texts: OpenLayers map interaction issue on margin left

texts:
I have created a sidebar that opens and closes. When it opens, it adds margin to the map container ( calling w3_open() ). See pen below.

However the interactions on the map e.g. clicking on the feature (getfeatureinfo) is wrong when there is margin on the left. I click on one feature but another feature is clicked

Probably I need to update map view when there it adds margin?

  Also note that when open and close browser's console for some weird reason the problem is fixed.

https://codepen.io/csandreas1/pen/xoWGyB?editors=1010

<script>
function w3_open() {
  document.getElementById("main").style.marginLeft = "25%";
  document.getElementById("mySidebar").style.width = "25%";
  document.getElementById("mySidebar").style.display = "block";
  document.getElementById("openNav").style.display = 'none';
}
function w3_close() {
  document.getElementById("main").style.marginLeft = "0%";
  document.getElementById("mySidebar").style.display = "none";
  document.getElementById("openNav").style.display = "inline-block";
}
</script>


texts: Batch uploading PDF files to ArcGIS Online?

texts:
I am trying to upload multiple pdf files to ArcGIS Online content all at once.

Is there is a script that I can use to add content to ArcGIS Online?


texts: Creating parallel line in ArcGIS Pro using ArcPy?

texts:
I want to create offset/parallel lines for my centre lines (over 12000 polylines) in different distance for each class number(i have numbers in column) and with the same line direction as centre line.

My idea was to create mid point for each line and then move the line at the same angle but i don't know how to do it in Python.

Also i know the way creating line buffer, splitting buffer etc. but then lines are not in the same direction and ArcGIS Pro basic not support polygon to lines tools etc.

The problem is that I am using basic licence in ArcGIS Pro which not support advanced tools.

Is there any way to do it using Python script?

import arcpy

#setup inputs
inFC = arcpy.GetParameterAsText(0)
outFC = arcpy.GetParameterAsText(2)

#setup output path
fcPath = outFC.rpartition("\\")[0]
fcName = outFC.rpartition("\\")[2]

if arcpy.Exists(outFC):
    arcpy.Delete_management(outFC)
arcpy.CreateFeatureclass_management(fcPath,fcName,"POLYLINE","#","DISABLED","DISABLED",inFC)
arcpy.AddField_management(outFC,"inFID","LONG","#","#","#","#","NULLABLE","NON_REQUIRED","#")

def shift_features(inFC):
    """
    Shifts features by an x and/or y value. The shift values are in
    the units of the inFC coordinate system.

    Parameters:
    inFC: string
        An existing feature class or feature layer.  If using a
        feature layer with a selection, only the selected features
        will be modified.

    x_shift: float
        The distance the x coordinates will be shifted.

    y_shift: float
        The distance the y coordinates will be shifted.
    """

    with arcpy.da.UpdateCursor(inFC, ['SHAPE@XY', 'Shift_X_right', 'Shift_Y_right']) as cursor:
        for row in cursor:
            cursor.updateRow([[row[0][0] + (row[1] or 0),
                               row[0][1] + (row[2] or 0)]])

    return

outFC = shift_features(inFC)

I try to shift my lines using attribute stored inside my table. When i put normal x_shift value like 2 meters shift function work correctly. But when i try to call my table it creates error. I don't know if shift function support shifting by attribute table values. Can you help me?

UPDATE
If that will help i got script that do the same us feature vertix to points.

# Feature Class to Points
#
# Paul Smith (2012) paul@neoncs.com.au

# Imports
import arcpy import numpy

#Inputs from user parameters
InFc  = arcpy.GetParameterAsText(0)
OutFc = arcpy.GetParameterAsText(1)

# Spatial reference of input feature class
SR = arcpy.Describe(InFc).spatialReference

# Create NumPy array from input feature class
array = arcpy.da.FeatureClassToNumPyArray(InFc,["SHAPE@XY"], spatial_reference=SR, explode_to_points=True)

# Check array and Exit if no features found if array.size == 0:
arcpy.AddError(InFc + " has no features.")

# Create a new points feature class else:
arcpy.da.NumPyArrayToFeatureClass(array, OutFc, ['SHAPE@XY'], SR)


texts:
I have created a model in Enterprise Architect (in ArcGIS technology). Inside I have several point, polygon, line feature classes and a couple of domains which correspond to them.

From the level of EA I export my model to ArcGIS XML Workspace. Then when I try to import it in ArcGIS 10 I get an error as below:


texts: Code only saving last polyline with ArcPy

texts:
My shapefile contains multiple polylines to segment and buffer and save in a new shapefile. My ArcPy code only saves the last polyline.
How do I save each output?
def main():
import arcpy

fc = r'D:\tayyab\New folder\buffer\cut_buffers\line.shp'
fc_out = r'D:\tayyab\New folder\buffer\cut_buffers\buffers_v11'
length_section = 100.0 # each 100 meters
buffer_dist = 10.0
tol_dist = 0.1 #
sr = arcpy.Describe(fc).spatialReference

# get polyline (first one from fc)
with arcpy.da.SearchCursor(fc, 'SHAPE@')as cursor:
    for row in cursor:
        polyline = row[0]
        lst_buffers = []
        start = 0
        end = length_section
        print "polyline.length", polyline.length
        while start < polyline.length:
            segment = polyline.segmentAlongLine(start, end, False)
            buffer1 = segment.buffer(buffer_dist)

            if end > polyline.length:
                end = polyline.length
            print " - start", start, "  end", end

            # get points around start and end of segment to determine angles
            pntg_d1, pntg_d2, pntg_start = getPoints(polyline, start, tol_dist)
            pntg_d3, pntg_d4, pntg_end = getPoints(polyline, end, tol_dist)
            pntg_center = polyline.positionAlongLine((start + end) / 2.0, False)

            # angle at start and end of segment
            angle_start = getAngle(pntg_d1, pntg_d2)
            angle_end = getAngle(pntg_d3, pntg_d4)

            # create perpendicular cutline
            cut_line_start = createPerpendicularCutLine(pntg_start, angle_start, buffer_dist, sr)
            cut_line_end = createPerpendicularCutLine(pntg_end, angle_end, buffer_dist, sr)

            # cut buffer polygon
            polygon = cutPolygonWithLines(buffer1, [cut_line_start, cut_line_end], pntg_center)
            lst_buffers.append(polygon)

            # for next segment
            start = end
            end += length_section

        # store results
        arcpy.CopyFeatures_management(lst_buffers, fc_out)

def cutPolygonWithLines(polygon, lst_lines, pntg_center):
    for line in lst_lines:
        try:
            lst_pols = polygon.cut(line)
        except:
            lst_pols = [polygon]
        polygon = getCorrectPolygon(lst_pols, pntg_center)
    return polygon

def getCorrectPolygon(lst_pols, pntg_center):
    correct_polygon = None
    for polygon in lst_pols:
        if polygon.contains(pntg_center):
            correct_polygon = polygon
            break
    return correct_polygon

def getPoints(polyline, dist, tolerance):
    d1 = dist - tolerance if dist - tolerance > 0 else 0
    d2 = dist + tolerance if dist + tolerance < polyline.length else polyline.length
    pntg_d1 = polyline.positionAlongLine(d1, False)
    pntg_d2 = polyline.positionAlongLine(d2, False)
    pntg_dist = polyline.positionAlongLine(dist, False)
    return pntg_d1, pntg_d2, pntg_dist

def getAngle(pntg1, pntg2):
    return pntg1.angleAndDistanceTo(pntg2, method='PLANAR')[0]

def createPerpendicularCutLine(pntg, angle, dist, sr):
    pntg_cut_1 = pntg.pointFromAngleAndDistance(angle - 90, dist * 2.0, 'PLANAR')
    pntg_cut_2 = pntg.pointFromAngleAndDistance(angle + 90, dist * 2.0, 'PLANAR')
    cut_line = arcpy.Polyline(arcpy.Array([pntg_cut_1.firstPoint, pntg_cut_2.firstPoint]), sr)
    return cut_line

if __name__ == '__main__':
    main()


texts: Calculating number of overlapping pixels across bands with non-zero values using QGIS

texts:
I have ~90 .tif files that are rasters of known fire boundaries from the last 90 years. I want to merge them so that each pixel has the value of the number of fires occurring at that pixel.
I can merge the raster files OK, but then each layer becomes a band, "nodata" sections become 0s, and the value of each pixel is the year (see table below from the identify function).

As you can see, this has "worked" in that I can click any pixel and see where there has been fires and when. What I now need to do is either sum or stylize this layer so that is shows the number of fire events (intervals not important). I cannot sum the value column as the value is the year rather than a "1" to correspond to a fire event.
Zonal statistics doesn't help me here because I need a viewable map rather than a table.
I have very little experience with raster-based GIS.


texts: Relating Tables in OLE DB from ArcGIS for Desktop?

texts:
I have a large registry of about 260,000 records  in microsoft access with data on land quarter sections. This data was to much for ArcGIS 10.2 to handle directly, so I created a Microsoft Jet 4.0 OLE DB connection to the database and added the query table I greated in Access.

Now I want to relate the table with the quarter section feature class (polygons of all the quarter sections). With this, I need to be able to select one or more quarter sections on the feature class and bring up all the related data from my access table.

However, ArcGIS help states that

  Because OLE DB connections do not have an ObjectID field, which is a
  field maintained by ArcGIS that guarantees a unique ID for each row in
  the table, you won't be able to perform certain operations in ArcMap.
  These include selecting the features in the layer on the map in any
  way, applying a definition query, creating full relates, or editing
  the attributes. See About ObjectID fields for more information about
  working with tables that lack Object IDs.

So I cannot relate a table from the OLE connection with ArcGIS due to lack of ObjectID.

Does anyone know a work around or a way to create the ObjectID so I can select related records from my large access dataset based on selected quarter sections on my feature class.


texts:
I have multidimensional meteodata I want to upload on GeoServer (2.17) using NetCDF format. So far, I don't care about projection other than WGS84. Following a documentation I tried to create in R test file with 1 variable and dimensions in the following order:
3 custom dimensions, lat, long (according to documentation I skip time and z or could I?). However, GeoServer logs showed that all 'custom' axes are non-recognizable
Unsupported axis: int here_is_name_of_my_custom_dim(here_is_name_of_my_custom_dim=9)

So general question - is it possible to skip time and z dimension and use only customs+2coordinate? If yes, then can anybody point to an example or post/tutorial/whatever to configure such a file (better with R, but Python also works). It's also ok if there is no way around to add 'empty' time and z dimensions (i.e. with 1 value), but it should work in the end.
P.S-  As an addition, I can't even reproduce example with time,z, lat, long dimensions (see R code below).
GeoServer gives error (copied only main part (as I guess, can parse all log))
java.lang.Integer cannot be cast to java.util.Date

library(ncdf4)
d1 = ncdim_def("time", "days since 2020-01-01 0:00:00", 1:5)
d2 = ncdim_def("z", "meters", 1)
d3 = ncdim_def("lat", "degrees_north", 1:10)
d4 = ncdim_def("lon", "degrees_east", 1:5)

v_def = ncvar_def('myvar', '', list(d1,d2,d3,d4), missval=-9999, prec='float')
nc = nc_create('test.nc', v_def, force_v4=T)
ncvar_put(nc=nc,varid='myvar',vals=1:250, start=c(1,1,1,1), count=c(-1,-1,-1,-1))


texts: Production Mapping "Load Data" Tool having issues with SDE

texts:
I am attempting to use the Production Mapping tool "Load Data" to import data from many personal databases into a new SDE database (SQL server RDBMS). The schema between the PGDB's and SDE is the exact same.

For some reason, I receive the error

  Target dataset [dataset name] not found or is inaccessible

although I am a dbo and can use the Data Loader or Append tool just fine to accomplish the same task (the Production Mapping tool allows multiple databases to be input at once, which is why I would love to get it working).

I correctly set up a cross-reference database between the databases and validated it, but still no luck.

Might this be a permissions issue? I am looking for a way to avoid tediously appending each dataset to the new SDE.


texts: Extracting data based on lat-long range using QGIS?

texts:
I have data on climate and was able to view it in QGIS. It shows me a map of the world with black/white color that's lighter or darker based on the factor my data was measuring. Since the data is for the whole world, it's really big and I can't open it on a data analysis site that was required of me.

Is there a way to limit the data to a certain lat-long range? (my data is lat-long coordinates with inf).

Could one search for NY or Mexico in QGIS and just get the data points for that region?

If this is not possible to make the data set smaller, would I just need to get data that's specifically for those places to start with?


texts: How to Populate new field using expression containing MAX value from another field in ArcGIS attribute table using field calculator?

texts:
I've got one field X with 10 numerical attributes (float/double).

I've created a new field Y (also float/double) which I need to populate with the percentage of each record/attribute (in each of the 10 rows) relative to the MAX value in field X, i.e. something like Y=X/max(!X!).

I need to do this using the field calculator - in VB script or Python.

Any suggestions? I thought using max before the field name would suffice? (I assumed it was a built-in-function from Python in ArcGIS).


texts: QGIS crashes when deleting memory layer created by a script

texts:
I wrote some scripts that create shapes based on different sources (Excel, XML...).
These scripts follow the template on Anita Graser's site and are integrated in my processing toolbox.

They work fine exploring, opening files and aggregating data as intended, and generating layers in memory.
Those layers can be saved and worked on.

But when I try to delete the generated memory layers, QGIS quits without an error message. I also can't use the "zoom to layer" function on these layers.

I may have forgotten something in my script but I can't find what.


texts: Labels and marker size in Print Layout don't match those on screen

texts:
I'm new to QGIS but have been flummoxed at almost the first turn.
My map looks great on the screen, but when I move to Print Layout and add a map, all the symbols and labels are much larger and spoil the map. Screengrabs:

All the sizes are, I think, in millimetres. I just use the defaults.
It's not practical to manually change sizes for every layer as I'm using the Natural Earth, which has lots of layers.
Is there a way to change symbol and label sizes globally? I've tried changing the paper size, but that doesn't help.
I'm using QGiS 3.10.10-A on a recent Macbook Pro.


texts: If I am using the coordinate system WGS84 ESPG:4326 and my cell size is 0.00833333,-0.00833333, how can I calculate the area of one cell in metres?

texts: I am currently trying to undertake some area calculations on some rasters. I have them all in the coordinate system wgs84 ESPG:4326 which I have just realised is in degrees. My cell size is 0.00833333,-0.00833333. How can I figure out the cell size in metres.


texts: Dividing space into homogeneous connected areas

texts:
Is there an algorithm (or software) for splitting a map into K connected areas such that the within-area variance of some variable is minimised, and the between-area variance of that variable is maximised?

This seems like it'd be a fairly common problem. For example, one might want to do this sort of thing if one wanted to draw homogeneous electoral districts.


texts:
How to open Cesium 3D Tiles with Qgis
The format is not available but this is an OGC standard:
https://www.ogc.org/standards/3DTiles

GeoPackage seems to be working on adding the extension:
https://www.compusult.com/html/OGC/3DTile_GeoPackage_Ext_Draft.html


texts: I need to add a sentinel-2 image into Google Earth Engine Explorer with a specific ID

texts:
I'm new at Google Earth Engine and Javascript.

I am using the graphical user interface of Google Earth Engine (Explorer) to obtain NDVI using sentinel 2. My study site is somewhere in Peru with an id:

S2A_MSIL1C_20180110T150701_N0206_R039_T18KYG_20180110T162736.SAFE

However, the format required is:

COPERNICUS/S2/20151128T002653_20151128T102149_T56MNN

and so I attempted to truncate my id using this code:

var image = ee.Image('COPERNICUS/S2/20180110T150701_20180110T162736_T18KYG');

But I get in response: Image.load:

Image asset 'COPERNICUS/S2/20180110T150701_20180110T162736_T18KYG' not found.

My complete code is:

//  image

var image = ee.Image('COPERNICUS/S2/20180110T150701_20180110T162736_T18KYG');

// Visualization parameters

var visParams = {bands: ['B8', 'B4', 'B3'], max: 3048, gamma: 1};
var visParams_ndvi = {min: -0.2, max: 0.8, palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
    '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'};

// Calculate NDVI

var image_ndvi = image.normalizedDifference(['B8','B4']);

// Map results

Map.centerObject(image,9)
Map.addLayer(image,visParams,'Sentinel-2 False Color Infrared')
Map.addLayer(image_ndvi,visParams_ndvi,'Sentinel-2 NDVI')


texts: Why does my "Zoom to selected rows" Shortcut not work in QGIS 2.6.0?

texts:
I´ve changed a shortcut, in this case "Zoom to selected rows". But now it stopped working. The new shortcut never worked and the originally as well is not working now. The icon gets activ when doing the shortcut, but no zoom-action. I´ve restarted QGIS, but without luck. Is this a bug?

WINDOWS 8 64-BIT
Standalone Install


texts: Can I automate using different reference grids in a data driven map

texts: I am trying to create a data driven layout that includes an overview map, at 1:15000, and a number of inset maps, at 1:4000.  I need to have two reference grids for the layout since the map scales are different and I'm trying to figure out if there is a way to switch between the two grids automatically.  Right now I am manually switching between the reference grids and then manually cycling through the pages to print them out.  Is this something that can be automated through the data driven pages UI or would I need to do some scripting?


texts: Using the output of zonal statistics in raster calculator

texts:
In ArcGIS (10.2) Model Builder, I am using the raster calculator to select a subset of cells from a temperature layer (a range between a max and min value). My expression reads as follows:
Con(("%temp%"  >= 7.25) & ("%temp%" <= 11.375), 1,0), where "%temp%" represents my temperature layer.
In plain English, I am asking it to select cells with a temperature value between 7.25 and 11.375.
The expression works fine, however, I will need to repeat this process multiple times, and each time the max and min temperature values will change.  Instead of inputting the values manually, I have been looking at calculating my temperature values with zonal statistics (this part works fine) and then using the output of the zonal statistics to give the max and min temperature values. For example, I have tried Con(("%temp%"  >= "%min%") & ("%temp%" <= "%max%"), 1,0), where "%min%" and "%max%" represent the output of the zonal statistics. Unfortunately, I am having no luck!


texts: How do I force the comma into Acres number on the ModelBuilder?

texts:
I'm writing a ModelBuilder to run many shapefiles output to database that has the sum acres for each crop however by looking at the field column it does not show comma  like for example

Alfalfa   351178.314628   what I would like to make the number 351178.314628 into 351,178 ?

I can export it to python and add if there is code; however I am not python expert.

I'll have this screenshot sample:


texts: How to arrange months properly in ArcGIS dashboard

texts: I am trying to edit my category selector for month in ArcGIS dashboard as it is displaying the months alphabetically instead of January, February, March and so on. How do I rectify this?


texts: Interrupting label display with ArcGIS Desktop?

texts:
In order to avoid slowness I would like to interrupt labels display when I'm working with ArcGis 10.

I read that the F9 key could do that but it doesn't work for me.
Does anyone know an other way to proceed, please ?


texts: ArcGIS .Net Framework SDK for 10.8: Does not install correctly

texts:
When using the installer:  ArcObjects_SDK_for_NET_Framework_108_172796
with having Visual Studio 2017 and ArcGIS 10.8 installed, the select features is missing the VS 2017 (see second image)
Steps already taken:

Uninstalled ArcGIS 10.8 and reinstalled it.
Updated VS 2017

Any suggestions?


texts: Python Error while installing QSpatiaLite in QGIS 2.0

texts:
I'd like to install QSpatiaLite to my QGIS2.0 but the following error massage occured while installing with the pluginmanager:

Couldn't load plugin QspatiaLite due an error when calling its classFactory() method

Traceback (most recent call last):
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 204, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "/home/map-pc/.qgis2/python/plugins/QspatiaLite/__init__.py", line 24, in classFactory
    from qspatialite import QspatiaLite
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 453, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/map-pc/.qgis2/python/plugins/QspatiaLite/qspatialite.py", line 29, in
    from QspatiaLiteApp import QspatiaLiteApp
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 453, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/map-pc/.qgis2/python/plugins/QspatiaLite/QspatiaLiteApp.py", line 27, in
    from Classes.Database import *
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 453, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/map-pc/.qgis2/python/plugins/QspatiaLite/Classes/Database.py", line 1, in
    from pyspatialite import dbapi2 as sqlite #Load PySpatiaLite
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 453, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/lib/python2.7/dist-packages/pyspatialite/dbapi2.py", line 27, in
    from pyspatialite._spatialite import *
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 453, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: libspatialite.so.5: Kann die Shared-Object-Datei nicht öffnen: Datei oder Verzeichnis nicht gefunden

Python version:
2.7.5+ (default, Sep 19 2013, 13:52:09)
[GCC 4.8.1]

QGIS version:
2.0.1-Dufour Dufour, exported

Python path: ['/home/map-pc/.qgis2/python/plugins/processing', '/usr/share/qgis/python', '/home/map-pc/.qgis2/python', '/home/map-pc/.qgis2/python/plugins', '/usr/share/qgis/python/plugins', '/usr/lib/python2.7', '/usr/lib/python2.7/plat-x86_64-linux-gnu', '/usr/lib/python2.7/lib-tk', '/usr/lib/python2.7/lib-old', '/usr/lib/python2.7/lib-dynload', '/usr/local/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages/PILcompat', '/usr/lib/python2.7/dist-packages/gtk-2.0', '/usr/lib/pymodules/python2.7', '/usr/lib/python2.7/dist-packages/ubuntu-sso-client', '/usr/lib/python2.7/dist-packages/wx-2.8-gtk2-unicode', '/home/map-pc/.qgis2/python/plugins/DigitizingTools/tools', '/usr/share/qgis/python/plugins/fTools/tools']


texts: Writing Where Clause for text field in arcpy.SelectLayerByAttribute_management()?

texts:
I am having problems writing WhereClause for my program.
I am using 10.3 arcmap.
My layer comes form a file geodatabase.

I have this code:

mlayer = arcpy.mapping.ListLayers(mxd, "1415_LGA_SpatJ", df)[0]
fn_field = 'MAPID'

....function which creates list of values a_list

for each_a in a_list:
    whereClause = '"%s" = %s'%(fn_field, each_a)
    arcpy.SelectLayerByAttribute_management(mlayer, "NEW_SELECTION", whereClause)

....

This code works when values of my field are numbers (field is long integral)

Now I need to modify the code and use another field instead, which has strings.
I get an error, invalid sql expression.

Could someone please help with updating sql statement.

I tried few things, but don't understand it really well.


texts: Adding a ruler to the OpenStreetMap site?

texts:
As a cyclist, I often need to measure the distance between two locations to have a rough idea.

The OSM website contains no ruler, and I was wondering if there were a way to add one, or if a similar site exists that does have one.


texts: PostGIS: Select FROM AS with WHERE not working

texts:
Maybe it's just because it's a Friday afternoon, but why can't I use the last WHERE in the following query? The ORDER BY works, but when I include the WHERE, I get an error saying "column line_geom does not exist".

SELECT ST_AsText(ST_Intersection(the_geom,ST_MakeEnvelope(...))) as line_geom
FROM (
  SELECT the_geom FROM table WHERE time > '2013-09-01'
) AS foo
-- WHERE GeometryType(line_geom) = 'LINESTRING'
ORDER BY line_geom

Thanks!


texts: Bleeding/fuzzy edges when adding raster to mosaic in ArcMap

texts:
I'm a relative newcomer to ArcMap.
I need to create one layer from three raster tiles that show presence/absence (0/1) of woody vegetation (background below). The three layers are publicly available here:
ftp://tern-auscover.science.uq.edu.au/nsw_spot_woody_extent_and_fpc/
The specific files are:
s5hgps_r419c105_y20082011_bcum5_r5m.img
s5hgps_r419c106_y20082011_bcum5_r5m.img
s5hgps_r420c105_y20082011_bcum5_r5m.img
Viewed by themselves, the raster layers are fine with sharp edges like this
.
After I create a mosaic dataset, and add the first raster layer it appears fuzzy on the edges like this.

I've been careful to match coordinate systems, and have tried snapping to the input file raster.
Any suggestions for what I'm doing wrong?
[the background: for my doctoral research, I need to measure the amount of woody vegetation within a 500m radius of my study sites on farms. My sites are spread across a large area that encompasses the three raster tiles above. Unfortunately, a few sites are directly on the join between two tiles (ie. they overlap). Therefore I am trying to create a single layer so that I can then perform the analysis of amount of woody vegetation cover]


texts: How to extracting values of 40 plots from different bands of raster image in R?

texts:
I have 40 plots(samples) and I want to extract corresponding values from different bands of RapidEye image. So at the end I want to have totally 200 values of 5 bands in the format of txt.
Till now I could open the image and different bands separately. I used:

>library(raster)
>library(rgdal)
>library(glcm)

>band1<- raster("D:\\MSc\\Imagedata\\Clipped\\Extent_bavaria.tif", band=1, values=TRUE)
>band2<- raster("D:\\MSc\\Imagedata\\Clipped\\Extent_bavaria.tif", band=2, values=TRUE)
>band3<- raster("D:\\MSc\\Imagedata\\Clipped\\Extent_bavaria.tif", band=3, values=TRUE)
>band4<- raster("D:\\MSc\\Imagedata\\Clipped\\Extent_bavaria.tif", band=4, values=TRUE)
>band5<- raster("D:\\MSc\\Imagedata\\Clipped\\Extent_bavaria.tif", band=5, values=TRUE)

>image_stack <- stack(band1,band2,band3,band4,band5)
>plot(image_stack)

I have my samples in the format of (.shp) and also (.xlsx) with coordinates. Now I want to find the value for each of my samples in each band of the image.


texts: I have GRASS6.5 SVN installed and wanted to go for SVN installation of Quantum GIS in Ubuntu 11.04. So as instructed in QGIS Wiki, I added the qgis repository to my sources.list file. Now I am able to see qgis1.7.0svn version in Synaptic but when I try to install it, GRASS is a dependency. How can I install only Quantum GIS and not GRASS from here.


texts: Using auto-completion/code-completion in OpenLayers 3?

texts:
I wonder how I can make use of code-completion with ol3?

In debugging-tools like firebug I can use the completion

but I haven't figured out how to use the completion in an IDE.
( for OpenLayers2 it was discussed here: https://stackoverflow.com/questions/17093764/enable-code-assist-for-openlayers-javascript-library )

Has anyone suggestions how to setup the OpenLayers3-code-completion / code assist in an IDE like AptanaStudio or Eclipse?