texts: Problem with QGIS raster interpolation from a point shapefile after duplicating the shapefile

texts:
Let me explain a problem I had with QGIS to see if somebody else has experience it.

I created a raster from a point shapefile through the Raster-> Interpolation menu with no problems. After that, I saved the point shapefile with a different name (same format), removed some points, saved and tried to create a new Raster the same way, but this time all values of the created raster were nans.

Puzzled, I saved a new copy of the original point shapefile, made no edits, checked that the metadata, attribute tables, etc. were exactly the same and tried to create a Raster, but again all nans were produced. My solution to the problem finally was to copy the shapefile (and all related files) through windows explorer, then add them to the project and finally create the Raster.

My question is, what was I doing wrong in the first case? What is the correct way of creating a copy of a shapefile? Can anyone replicate the problem?

I am using QGIS 2.6.0 Brighton 32 bit for windows.

Regards.


texts: How to select highest value from multiple attributes' columns values in ArcGIS Model Builder?

texts:
I am using ArcGIS version 10.3. I have two shapefiles (points and polylines) and I am trying to create a model in ArcGIS Model Builder to automate my task. What I have done so far is:

Created new geometry around points using a buffer of 30m
Intersected the new geometry with lines layer to select lines (streets) within 30m buffer
Used an iterator to repeat this task for all points (87 points total)

Here is the screen shot of the model:

and here is output of the model.

and screen shot of output's attributes.

So far, I can view the address locations and buffered streets like this.

I have explained the layers as well. The scenario contains address and streets within 30m location. Each line feature contains the values of Annual Daily Traffic (ADT) from 1960 - 2005 (10 columns). The image above only shows the different ADT values of each line feature (using query tool). Can someone help me to?

Select the nearest street with highest ADT value (there are 10 ADT columns in each line feature)?
Shift my address location to that layer containing highest ADT value (this should create another layer for all addresses)?

Can someone help me to do that?


texts: When in editing session and merging polygons, feature labels for each are not from desired field in ArcMap

texts:
In my mxd, the polygons in the layer that I'm editing display a label from a certain chosen field. But when I am editing and merging these polygons, the label displayed in the 'Merge' dialogue box is from a different, less relevant field.
How do I make the 'Merge' dialogue box display the labels as shown in that polygon layer?


texts: Problem with MCC-LiDAR (insufficient memory)

texts:
I am running mcc-lidar on a win7 64bits with 16gb of ram. Nevertheless when I try to process a las file with 1, 800 000 points using -s 0.5 and -t 0.3 the program tells me:

Error: insufficient memory; if applicable use a larger post spacing.

I am able to process the data using -s 1.5 -t 0.3, so I figure the problem is not with the number of points.

The program is installed at c:\program files (x86), which I understand is for 32 bits programs. Is there some specific instructions to install it as a 64 bit application in order to use my 16 gb of RAM?


texts: PostgreSQL JDBC prepared statement query string with single quotes giving error

texts:
I'm trying to run a query using JDBC but i'm having difficulty injecting values into the prepared statement. Here is a sample of what I was doing:

String queryString = "... WHERE location <@ box '((?, ?),(?, ?))' ..."
PreparedStatement ps = this.connection.prepareStatement(queryString);
ps.setDouble(1, x1);
ps.setDouble(2, y1);
ps.setDouble(3, x2);
ps.setDouble(4, y2);
ps.executeUpdate();

Which gives me this error:

org.postgresql.util.PSQLException: The column index is out of range: 1, number of columns: 0

I think it thinks that the values in the single quotes are string literals and so it doesn't see the ? as parameters to inject. I've tried various other

Does anyone know how I could fix this? Or rather what else I should be looking to do?


texts: Invert legend values for raster in QGIS 3.4.6

texts:
The legend in the TOC for a raster layer ( > 1 million cells with continuous data) in my view might look like this:

The values go from low (above) to high (down). Is there a way to show the values in reversed order: high values on top and low values below?


texts: Creating buffer around points belonging to object using lidR

texts:
I would like to create a 3D buffer around powerlines inside a LAS file that contains forest and power lines in lidR (R package). If I would manually delineate the line features and export to an individual file,
How do I create a buffer around all points in that file?
I am thinking of something like this:


texts: GRASS - Calculate 90th percentile of a timeseries dataset

texts:
How can I calculate the 90th percentile from a registered timeseries dataset of 30 rasters in GRASS?
GRASS timeseries has 2 "statistical" commands which are t.info and t.rast.univar. Unfortunately, the first one calculate, from the whole set, which is what I want, only the maximum and minimum values, not percentiles/quantile.
On the opposite side, t.rast.univar calculates the percentile, togheter with other statistics (max, min, average and others) but what it does is to calculate it for every map registered, not for the whole set.
Is there a way to "integrate" both commands and to calculate the 90th percentile within all raster maps?
There is a third command outside the timeseries package to analyze statistics, r.series, but it calculates quantiles only for the single cell within all input maps. What I want is the 90th percentile of the whole area of the 30 rasters I have (which is global).


texts: ArcGIS Pro 1 year home use says I'm not authorized

texts:
I just downloaded and installed the ArcGIS 10.3 and ArcGIS Pro 1 year home use. ArcGIS 10.3 is working fine. However, I have to sign into esri online when I initiate ArcGIS Pro. When I sign in using my esri account I bought the product with, I get this:

I can't seem to locate any information in my install email or online about how to register for ArcGIS Pro. I'd think that since it's bought and paid for through the account it would be connected.

Am I missing something?


texts: Get latitude and longitude of nearest seaside

texts:
I have latitude and longitude of a point and I want to find out lat&long of (or distance to) a nearest seaside point.

Which way to do this will you recommend?

GIS technologies are new to me.


texts: Split vector layer into GeoJSON files as a batch

texts:
I have a large set of polygons in shp, each with unique string value in attribute table. Now I need each of them in separate .geojson file (one polygon = one .geojson) but with name equal to the value in table.

I've found Multipart to singleparts in default processing panel which works just fine with the iterate option enabled, but there's no way to choose name from table. I have never wrote a script for qgis but I have some basic knowledge of programming and editing code so I tried to find this script on the computer and modify it. However I can't locate the file by name.

Or maybe there is another option to write script using this ready piece? Or some plugin which allows saving to multiple files with chosen filename?

I'm using QGIS 2.18 and QGIS 3.0.


texts: Getting additional subbasin areas designated by points in river network？

texts:
I follow Watershed and Stream Network Delineation, and here's my workflow：

use Batch Point Generation to create batch points related to designated points(green dots)
Batch Subwatershed Delineation using these batch points

and I got this：

It looks great, but I still can't clip existing subbasins using designated points to get smaller subbasins inside.

Are there tools in Arc Hydro Tools or any other tool like QGIS or GRASS that can help me get my goal？

I have got the subbasin area of every river segment using Arc Hydro Tools：

Martin suggested How to get the subbasin area of each specified river segment？

Now I have some sample stations(the green dots) in the river network, and I'd like to take them as the segment bases, thus I would get more segments, and I want to split the drainage line based on those points and calculate the subbasin area of these new segments：

I think the key is how to make ArcGIS take these green dots as confluences of rivers as well. And all I want is like this schematic photo：

What can I do?


texts: Does Collada allow to georeference the 3D coordinates within a DAE file?

texts: Can I assign lat/long/height values to each vertice of a Collada 3D model?


texts: Are there limitations using a PostGIS out-db raster?

texts:
I have been trying to import a large raster into a PostGIS 2.0 database and hit a few stumbling blocks.

Initially when I tried to import the raster using raster2pgsql, I would receive the following error:
rt_band_set_pixel_line: Coordinates out of range

After finding this page, and using gdalinfo I found that the raster exceeded the maximum allowed size of 65535 x 65535.

gdalinfo output:

  $ gdalinfo ari100.tif
  Driver: GTiff/GeoTIFF
  ...
  Size is 42971, 77138
  Coordinate System is:
  GEOGCS["WGS 84",
  ...
      AUTHORITY["EPSG","4326"]]
  Origin = (152.741676637167842,-26.215302802012008)
  Pixel Size = (0.000009239757419,-0.000009239757419)
  Metadata:
    AREA_OR_POINT=Area
  Image Structure Metadata:
    INTERLEAVE=BAND
  Corner Coordinates:
  Upper Left  ( 152.7416766, -26.2153028) (152d44'30.04"E, 26d12'55.09"S)
  Lower Left  ( 152.7416766, -26.9280392) (152d44'30.04"E, 26d55'40.94"S)
  Upper Right ( 153.1387183, -26.2153028) (153d 8'19.39"E, 26d12'55.09"S)
  Lower Right ( 153.1387183, -26.9280392) (153d 8'19.39"E, 26d55'40.94"S)
  Center      ( 152.9401974, -26.5716710) (152d56'24.71"E, 26d34'18.02"S)
  Band 1 Block=42971x1 Type=Float32, ColorInterp=Gray
  ...

So I decided to use the out-db functionality by using the -R flag in raster2pgsql.

From the documentation:
Register the raster as a filesystem (out-db) raster, Only the metadata of the raster and path location to the raster is stored in the database (not the pixels).

I imported the raster using the following command:

raster2pgsql /gis/Flood/ari100.tif -R | psql -U username database

But when I tried a query:

SELECT ST_Value(rast, ST_PointFromText('POINT(152.9632 -26.4878)')) FROM ari100

I get the following error:

NOTICE:  Attempting to get pixel value with out of range raster coordinates: (23975, 29491)
CONTEXT:  PL/pgSQL function "st_value" line 13 at RETURN

However, if I use GDAL (through Python) to query the raster directly using this script

python val_at_coord.py 152.9632 -26.4878

I can retrieve the value correctly.

I decided to find the extents of the raster in the database:

SELECT ST_Height(rast) As rastheight, ST_Width(rast) As rastwidth from ari100;

rastheight | rastwidth
11602      | 42971

The height is completely incorrect (should be 77138).

So my question is: did I do something wrong (in using raster2pgsql) or is there a limitation in using out-db rasters with PostGIS?

Thanks


texts: How to allow more polyline vertices for a ST_GEOMETRY in Oracle/SDE?

texts:
I'm receiving an error writing geometry (in ArcObjects) that has a large number of vertices (around 900).  So big, but not ridiculously big.  I would expect this to be no problem and I don't recall this issue in SQL Server or PostGreSQL.  The geometry storage is ST_GEOMETRY.

The error is:
ORA-24816: Expanded non LONG bind data supplied after actual LONG or LOB column

..and it appears to be related to the number of vertices that the polyline has on it.  Lots of vertices: I get the above-mentioned error.  Less vertices: no problems.

As a test, I generalized the geometry down to two points: success.

I've started to work my way up and gotten up to 120 vertices with success.  Is there a better way to determine that limit?  Is there a way to increase the number of vertices that I can have on a line?  Is there a DBTUNE parameter that manages this?   I've googled the actual ORA error, but I'm uncertain as to how to apply this to an ArcSDE repository.

A couple of notes:

I don't own or administer the ArcSDE repository and can't modify it.
Breaking the long polyline feature into several smaller ones is not
an option (for this project).
Write using feature.Shape = geometry; feature.Store();
Code executes as expected for most of the features--I have 2 that have high number of points that cause the ORA error.

Oracle version:

Oracle Database 11g Enterprise Edition Release 11.2.0.2.0 - 64bit Production
PL/SQL Release 11.2.0.2.0 - Production
CORE    11.2.0.2.0  Production
TNS for Linux: Version 11.2.0.2.0 - Production
NLSRTL Version 11.2.0.2.0 - Production


texts: Comparing different open source GIS servers?

texts:
I currently know the following open source GIS servers:

Geoserver: java based, vector and bitmap support, also supports. Is now on par with mapserver (performance-wise).
mapserver: c/c++, used to be the fastest server. Better for raster than vector?
mapguide: do not know this
mapnik: the new kid on the block? do not know much about it, but it looks appealing
mapfish: as far as i know, only retrieves vector-data for display in openlayers. There does exist a rails implementation.

If possible i would like some kind of comparison, which did you choose or prefer and why?

I am looking at building a rails website, and need some kind of GIS server. I will need raster and vector data (clickable). Is this doable with one server? Backend database will be Postgis.


texts: QGIS 2.18 - keeping field properties when using data management or copying layers

texts:
I want to use Data management tool to "join attributes by location".  this creates a temporary layer with the data.
I can apply the style of the original layer, but not the field properties.  Is there a way I can either add the new joined information to the original layer, OR apply the field properties to the temporary layer when I save as?
this is version 2.18, aware it is an option in later versions


texts: Join CSV file to shapefile using gdal/ogr

texts:
I have a shapefile with several attributes, for example "YEAR", "COUNTY", and "AREA". I also have a CSV file with more fields that I want in the shapefile, such as "POPULATION". Both the shapefile and the CSV file have a field "GISJOIN".
I know how to do a join in QGIS. But how can I make a permanent join and write to a shapefile using ogr2ogr or one of the other tools in GDAL/OGR?


texts: I have installed PostGIS 2.0 and need new topology but it doesn't seem installed by default. I found a lot of command line to install PostGIS 2.0 with topology but is there a way just to add topology?


texts: calculate distance between points in different point layers

texts:
I have two point layers.
Layer one contains a single point (a landmark)
Layer 2 contains a set of points (person locations)

Is there an automated means of obtaining the distance of each individual point in layer 2 from the single point in layer 1.  And can this result be filled to a field in attribute table of layer 2.

(Complete novice!)


texts: How to enter formula into ArcGIS Raster Calculator?

texts:
i have the following formula:

L = (l/22,1)^m

(l = Flow length)

m is calculated like this:

m = 0,5 when slope is >4 degrees

m = 1,2 * sin^1/3a  when slope is <4 degrees   (a = slope in degrees)

I already have the slope_deg and flow_length layer. Now I want to enter the formula to the raster calculator. I am new into ArcGIS and I dont know how i can implement it.


texts: Using gdal.Warp() in Python to clip a raster to a shapefile is producing errors

texts:
I have some super simple code that I'm using to clip raster data to US Highways. However, the clip isn't producing pixels for everywhere that the shapefile intersects the raster and I cannot figure out why this is the case. Here is the code:
gdal.Warp(input_raster_filepath, output_raster_filepath, cutlineDSName=road_buffered_filepath, cropToCutline=True)

Here is the raster and shapefile before the clip:

And here is the raster after the clip:

I was expecting for there to exist a pixel from the raster at every point along the highway but this is not the case. Can anyone explain why? Or could anyonen provide a better method for doing this?


texts: How do you pull EPSG:4326 tiles from GeopackageJS using OpenLayers 6?

texts:
I'm trying to leverage GeopackageJS for use in a mapping application using OpenLayers 6, but the documentation seems to be really iffy-- the link to the API documentation actually 404s and the examples on the actual github seem to be woefully out of date. I've had to build most of my understanding using the geopackage.js script itself as well as some helpful Stack Exchange answers.
The map I'm displaying is in an EPSG:4326 projection and I have the current solution using one that I found from this GIS Stack Exchange answer. I have a custom tile loader that essentially looks like this:
var promise = geopackage.getTileFromXYZ(this.gpkg, this.gpkg.getTileTables()[0], x, y, z, 256, 256);
promise.then(function(base64) {
    tile.getImage().src = base64
});

The problem is that getTileFromXYZ specifically has the TileRetreiver fetch the tile in EPSG:3857. So either the tile gets stretched and distorted to fit onto the EPSG:4326 map I have or it doesn't line up with the projection. The release notes for GeoPackageJS explicitly say they added a method to retrieve tiles in EPSG:4326 in version 1.1.4 but I cannot find any method that actually exposes that kind of logic.
Does anyone know what I'm missing or if there's a good example page or documentation I should be looking at?


texts: I am a newbie, I think I'm pretty resourceful, but I cannot figure out, or find out by searches, just what "map units" are in QGIS as compared to "decimal degrees". For example I can only navigate by entering coordinates in the status bar with "map units" - decimal degrees do not seem to work? I assume I am missing something? And is there a tool to convert "decimal coordinates" into "map units?"


texts: ArcGis Runtime SDK for .Net - keep feature in the middle of the screen

texts:
i'm working with ArcGis Runtime 10.2.5 SDK for .Net.

I want to keep a graphic in the middle of the screen no matter what happens in basemap. I want to simulate an aircraft flying over the map, keeping it pinned at the center and updating basemap envelope.

My first approach was update both basemap envelope and aircraft geometry with same the position, but since the procedures don't execute at same time

    //first basemap is centered
    MapPoint center = new MapPoint(aircraft.longitude, aircraft.latitude, SpatialReferences.Wgs84);
    MapView.SetView(center);
    //after that, aircraft is centered
    draw_aircraft(center.X, center.Y);

there are some "vibrations" when zoom is near to earth (where low difference in geo coordinates means big differences in pixels).

is there any way to do this?


texts:
I'm using an undocumented private API to grab some tiles. But it's using a strange coordinate system.

The area covered is

-77.318467, 44.148949, -77.297040, 44.133303

however the APIs bounding box is

-8607208.443178354, 5485854.929815598, -8604638.248102197, 5489027.066489516

Is there a name for this system? And is there a simple way to convert lat/lng to this coordinate system?


texts: Clipping and reprojecting in one process for modelbuilder // "Save raster layer as" ... in modelbuilder

texts:
I want to build a model in the QGIS Modelbuilder.
First it should do a clip. Afterwards it should be reprojected (from Geographic WGS84 Coordinates to Gauss-Krüger). I only see one way to do this: with the clip by extent and warp gdal tools. This will produce lots of missing values at the edges, which is quite bad for my application. Although it would work with another clipping, it is kind of complicated to ask for 8 coordinates in 2 different coordinate systems (which couldn't even be simply translated) and kinda circuitous for an automated process.

The "save as" / "save raster layer as" menu (by rightclicking the layer) would be optimal for this task, because it can directly extract in a (re-)projected Coordinate System. Is there any way to use it in the modelbuilder? I guess only by implementing it as an extra module. I can't find its source code, though. Do you know where I can find it?

Did I miss something / is there a better way to do this task?


texts: I am attempting to use the Grouping Analysis tool in the Spatial Statistics (in ArcGIS 10.2)  to group state high schools (points) together based on enrollment and location.  Is there a way to make each cluster the same size?  I am able to set how many groups I want, but I cannot seem to make each the same size.


texts:
I am running ArcMap 10.6 and I need a line that I thought would be included in the standard styles. It is the line symbol used when a polygon feature crosses a line feature and you want to show that on a map without labelling both sides.

In the form I have usually seen it, it is simply a thin line with short lines on opposite sides of both ends at acute angles to the main line.

Do any of you have that in your style sets?

Illustration of what I want

I want to create this using the drawing tool. I would probably then convert it to a feature.

How difficult is it to create new symbols for repeated use?


texts:
I've upgraded to GeoServer 2.3.0 running on Apache Tomcat 7.0.39 on my host WebFaction. My websites and GeoServer are on same host but different ports. How do I go about setting a proxy so that my site domain can grab data from GeoServer. I remember using a proxy.cgi with openlayers, but I'm now using Leaflet.

Is it possible to have my sites with domain capecodgis.com retrieve data from tomcat.capecodgis.com/geoserver.

I apologize in advance if I'm not clear. This is very confusing. I've been reading different docs and forums to no avail


texts: Count Overlapping Lines using ArcGIS Desktop?

texts:
I'm trying to calculate the frequency of overlapping lines.

The question is pretty identical to the illustration show here Frequency of partially or fully overlapping lines

However I have no python programming skills, so any advice or guidance on achieving a similar result would be a great help.

I'm using ArcGIS 10 Basic.


texts: QGIS Plugins - How to make sure that Vector and Raster layers in a folder are loaded only once?

texts:
Plugin Intro:

I am creating a plugin for in-house development which lists the location of vector/raster layers of specific projects. Each project has a number of files ranging from 0 to few hundreds. The preface of the plugin is, I need to type the project ID and click a button which fetches the list of files (from an in-house database) in the form of table widget. When I click a specific row and click another button the files from the row clicked are opened in QGIS.

Problem:

When the number of files is one, the plugin works perfectly. The problem arises when the number of files are more than one. In the latter case when I click a row the same file is loaded several times. There is no loop in my code so I don't understand why this is happening. The code is pasted below.

Function to retrieve data and show in table widget:

def ListPage(self):
    self.tableWidget.clear()
    self.cid = self.textEdit.toPlainText() #to get the number of project id #
    self.data = prj.getProductsByStatusAndOutput(status='READY_TO_CLOUDMASK', output_id=self.cid)
    self.tableWidget.setHorizontalHeaderLabels (('Id','File Name','Status', 'Location'))
    header = self.tableWidget.horizontalHeader()
    header.setSectionResizeMode(0, QtWidgets.QHeaderView.ResizeToContents)
    header.setSectionResizeMode(1, QtWidgets.QHeaderView.ResizeToContents)
    header.setSectionResizeMode(2, QtWidgets.QHeaderView.ResizeToContents)
    header.setSectionResizeMode(3, QtWidgets.QHeaderView.ResizeToContents)

    for rowIndex, record in enumerate(self.data):
        itemId = QTableWidgetItem(record['identifier'])
        itemStatus = QTableWidgetItem(record['status'])
        itemFilename = QTableWidgetItem(record['filename'])
        itemloc = QTableWidgetItem(record['prod_location'])

        self.tableWidget.setItem( rowIndex, 0, itemId)
        self.tableWidget.setItem( rowIndex, 1, itemFilename)
        self.tableWidget.setItem( rowIndex, 2, itemStatus)
        self.tableWidget.setItem( rowIndex, 3, itemloc)

        self.pushButton_2.clicked.connect(self.OpenQgis)

Function to open the files when a row is clicked:

def OpenQgis(self):

    row = self.tableWidget.currentRow()
    path = (self.tableWidget.item(row,3).text())
    sname = (self.tableWidget.item(row,1).text())
    vlocation = path + "/" + sname
    tname = (os.path.basename(vlocation)[0:28]+"_stack.tif")
    tfname = (os.path.basename(vlocation)[0:28]+"_stack")
    rlocation = path+ "/" +tname
    self.textEdit_3.setText(tname)
    if os.path.exists(vlocation):
        iface.addVectorLayer (vlocation,"","ogr")
    else:
        QMessageBox.about(self, "Error", "Invalid shape file location")
    if os.path.exists(rlocation):
        iface.addRasterLayer (rlocation,tfname,"gdal")
    else:
        QMessageBox.about(self, "Error", "Invalid raster file location")


texts: How to sum values of points inside polygons?

texts:
I have a polygon layer of territory and point layer with values from measurements.
I must to add sum values from points layer to polygon layer .
How to do this in PostGIS?

Supplemented:

Point layer is:

id|name   |value
1 |first  |547051
2 |second |539949
3 |third  |539864
4 |fourth |539724
5 |fifth  |500208

Pattern layer:

and result must be

id|name|sum
55|I5  |1619037
61|H7  |1047259


texts: Earth Engine authentication issue: "Request had invalid authentication credentials"

texts:
I am trying to load one surface reflectance image from Sentinel-2 MSI Level-2A in Google Earth Engine, but I get the following error:

  ImageCollection (Error)
  Request had invalid authentication credentials. Expected OAuth 2 access token, login cookie or other valid authentication credential. See https://developers.google.com/identity/sign-in/web/devconsole-project.

Here is my code and the roi:

var roi = ee.Geometry.Point(-94.86,36.93);
var image = ee.ImageCollection('COPERNICUS/S2_SR')
    .filterBounds(roi)
    .filterDate('2017-05-01', '2017-06-30')
    .sort('CLOUD_COVER')
    .first();
Map.addLayer(image, {bands: ['B4', 'B3', 'B2'],min:0, max: 0.2}, 'True colour image');

Any idea why this is?


texts:
I'm trying to get a list of my datasets from the API.  In version 1 of the API, a 'viz' call would return vizualizations (marked as "derived") and the underlying data (marked as type "table"), but it appears the API 1.0 isn't supported anymore, and version 2 doesn't return it.

I had hoped that the map call would returned a reference to the underlying data set(s), but that doesn't seem to be.

Any suggestions on how to programmatically get a list of an account's datasets?


texts: Data from google maps or bing maps available for offline use? Whole cities, I mean

texts:
I was developing some graph-based search algorithms and it would be interesting to use them against some real data. I used openstreetmap against them (downloaded some city-wide OML files from some obscure site), and it all worked well until I started noticing that the data wasn't that trustworthy -- the data near my house had wrong highways directions, and there were roads that weren't properly connected.

I was wondering whether I could download full-city google or bing (or actually, any found to be trustworthy source) maps data for this same batch processing? What I basically need is access to all the streets (and to know whether they are one or two-way streets) of a given town of my choice (Lisbon, for instance), as a combination of vertex + edges, or any other format that will let me gather that info, for offline use.


texts:
I just changed the mountain peak label from the outdoors style (mapbox).
As the default there is just a mountain icon + the name of the mountain.
I also wanted to label the elevation but when I import the style to QGIS my style change is not getting displayed.

Here the label when I import the changed style to QGIS:

Here the changes I made (displayed in the mapbox viewer) the label should look like this in QGIS:


texts: Create basemap tiles for a custom shape (for reducing the tiles count)

texts:
in order to create a world wide tiles (for a basemap) in acceaptable zoom level one must have a huge hd (see in this post). but, most of the world is sea...

so i want to create x levels of the whole world (x is small say 5-10) and then y levels of only the mainland (y is big say 20) and by that save allot of cpu time for the tile creatiion and hd space for saving unneeded tiles.

i use tilemill to create tiles and you can only define a bounding box rectangle.

can anyone recommend a solution (or maybe a different approach?)

thanks!

edit:
the tiles are created from osm data


texts: Why does Raster Calculator (RasterCalculator_sa) give ERROR 000539 in ArcPy?

texts:
I am getting a not executed error for a raster calculator.   Here is part of the output:

  File "C:\projectTemp\FloodAnalysis\floodxsectionanalysis.py", line
  107, in makeSeemless
           arcpy.gp.RasterCalculator_sa(Con(IsNull(rasterIn), FocalStatistics(rasterIn,NbrRectangle(3, 3), "MEAN"), rasterIn),
  "{0}yrSeemless.img".format(increment))   File "C:\Program Files
  (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\geoprocessing_base.py", line
  498, in
      return lambda *args: val(*gp_fixargs(args, True)) arcgisscripting.ExecuteError: ERROR 000539: Parsing error SyntaxError:
  invalid syntax (line 5) Failed to execute (RasterCalculator).

This code I have setup to put into this calculator is:

def makeSeemless(outputPath, increment, rasterIn):
    arcpy.gp.RasterCalculator_sa(Con(IsNull(rasterIn), FocalStatistics(rasterIn,NbrRectangle(3, 3), "MEAN"), rasterIn), "{0}yrSeemless.img".format(increment))


texts: Use QGIS or QField with a GPS physically walk a polygon in the field

texts:
Is it possible to use QGIS or QField with a GPS to physically walk a polygon in the field?

I'd like to use the software to physically walk shapes in the field but have been unable to find information on how to do so.


texts: ExecutING SQL script in ArcGIS field calculator?

texts:
I am trying to transfer over row data from a table to columns in a shapefile table. There is a common "SA" variable between the two (i.e. the study area number) but there is more than one land use type, hence why I must create two new columns for each landuse, in order to capture the area and percentage value for each type. A join won't work because it only captures the values for one land use as opposed to all of them, as the first table doesn't have as many entries as the second.

The reason I need to migrate the data is because I'm trying to use this in dynamic text on data driven pages.

I wrote a SQL script with the help of my friend but it won't execute through field calculator, giving the generic 99999 error:

UPDATE SApeelfinal
SET LUp1 = (SELECT perc FROM stat_fin WHERE stat_fin.SA = SApeelfinal.FID_text and LANDUSE = 1)
WHERE EXISTS (SELECT 1 FROM stat_fin WHERE stat_fin.SA = SApeelfinal.FID_text and LANDUSE = 1)

Though he said that I have to change the LANDUSE value for each time I run the script and as well change SUM_AREA and perc. I tried putting it in a toolbox but it wouldn't run that way either. Can someone enlighten me as to a better way to do this?

Also I was wondering if I had to manually create the new fields for this to work.

Here is a screenshot of the two tables:

You can also find a bigger version here: http://i.imgur.com/sx08u.jpg

The red is the common variable between the two tables, the blue is the field I want the percentage to go to (and there would obviously be a new one for area), and the green are the two fields I want to migrate over based on the SATXT and FID_text common fields.

I've tried using the following python script in the codeblock:

    def getNewVal(newVal,landuse):
  if landuse == 1:
   return newVal
  else:
   return

With this in the calculate field:

getNewVal(!stat_fin.perc!,!stat_fin.LANDUSE!)

But it returns the same error. Is there something I'm doing wrong?

To clarify, there are no joins because joins make dynamic text not work.


texts: GEE - Image collection to multiband image

texts:
I'm new to GEE and got stucked running the followed code. My problem is: Get L5 collection, selecting B7 channel > Mask clouds and shadows based on pixel_qa product > Convert the image collection to a multiband image (only B7) > Export to drive.

The code is running but I can't display the layer or export it. Error message: Cannot add an object of type ComputedObject to the map.

Does anyone know how can I do this?

// Test area
var geometry = /* color: #98ff00 */ee.Geometry.Polygon(
    [[[-49.7, -19.3],
      [-49.7, -18.9],
      [-50.1, -18.9],
      [-50.1, -19.3]]]);

// Function to mask clouds
var maskClouds = function(image){
  var pixel_qa = image.select('pixel_qa');
  return image.updateMask(pixel_qa.lt(69))
              .clip(geometry);
};

// Import L5 collections
var collection5 = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
    .filterDate('2008-01-01', '2008-12-31')
    .filter(ee.Filter.eq('WRS_PATH', 222))
    .filter(ee.Filter.eq('WRS_ROW', 73))
    .select('B7', 'pixel_qa')
    .filterBounds(geometry).map(maskClouds);

// Image collection to Image
var empty = ee.Image().select();
var landsat5 = collection5.iterate(function(image, result) {
  return ee.Image(result).addBands(image.select(['B7']));
    }, empty);

// Display
print(landsat5)
Map.addLayer(landsat5, {min:0, max:1000})
Map.setCenter(-49.9,-19.1, 10);

// Export image
Export.image.toDrive({
  image: landsat5,
  description: 'landsat5',
  region: geometry,
  folder: 'test',
  maxPixels: 1e12,
  scale: 30
});


texts: How to automate generating Directory of spatial files (shapefiles) and publishing contents via REST in GeoServer 2.14

texts:
I have a folder containing shape files for a few dozen layers.
On my development server, I created a Directory of Spatial files store in my workspace, then published each one, one at a time, with the UI.  Create the layer, select the store, click publish next to the next shapefile name, click calculate bounds, click calculate lat/long bounds, click Publish tab, set default style, Submit. Repeat for each layer.

Now I'm setting up a production environment.  I've installed GeoServer, and would like to import all these layers.  I've never used the REST API's, but am quickly coming up to speed.  I've used CURL / REST to create my workspace.  Now I'm trying to create my datastore.  I haven't found an example to setup a "Directory of spatial files" datastore.  I queried my existing store, so I can see the XML (or JSON) to post, but I get nothing as a response (I use -o option of CURL, but the output file is empty).

Here is the CURL command I use:

 curl -u admin:geoserver -v -XPOST -H 'Content-type:text/xml' -T My_datastore.xml -o CreateMyDataStore.out http://localhost:8080/geoserver/rest/workspaces/MyWorkspace/datastores

The My_datastore.xml file looks like this:

<dataStore>
    <name>MyDatastore</name>
    <type>Directory of spatial files (shapefiles)</type>
    <connectionParameters>
        <entry key="charset">ISO-8859-1</entry>
        <entry key="filetype">shapefile</entry>
        <entry key="create spatial index">true</entry>
        <entry key="memory mapped buffer">false</entry>
        <entry key="timezone">America/New_York</entry>
        <entry key="enable spatial index">true</entry>
        <entry key="namespace">http://www.opengeospatial.net/cite</entry>
        <entry key="cache and reuse memory maps">true</entry>
        <entry key="url">file://C:\Geo Server\MyFolder</entry>
        <entry key="fstype">shape</entry>
    </connectionParameters>
</dataStore>

The file directory "C:\Geo Server\MyFolder" is the same name used to create the datastore on the development machine...

The output I get at the command prompt when I run the CURL command is as follows:

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0*
  Trying ::1...
* TCP_NODELAY set
* Connected to localhost (::1) port 8080 (#0)
* Server auth using Basic with user 'admin'
> POST /geoserver/rest/workspaces/MyWorkspace/datastores HTTP/1.1
> Host: localhost:8080
> Authorization: Basic YWRtaW46Z2Vvc2VydmVy
> User-Agent: curl/7.61.1
> Accept: */*
> 'Content-type:application/xml'
> Content-Length: 614
> Expect: 100-continue
>
< HTTP/1.1 100
} [614 bytes data]
* We are completely uploaded and fine
< HTTP/1.1 415
< X-Frame-Options: SAMEORIGIN
< Access-Control-Allow-Origin: *
< Accept: application/json, text/json, application/xml, text/xml
< Content-Length: 0
< Date: Thu, 25 Oct 2018 03:25:30 GMT
< Connection: close
<
100   614    0     0  100   614      0   1710 --:--:-- --:--:-- --:--:--  1710
* Closing connection 0

I see the Error 415, but it doesn't make sense.  I've specified "text/xml" and it is one of the accepted content types.

Anybody see anything subtle or glaringly obvious that I missed?
Can I not create this type of store through the REST API?
I'd love to be able to automate this so I can fully automate the install/configuration of new production GeoServers (as they run locally at the customer site).


texts: Where to find street names for each city in US?

texts:
I'm sorry if my question is stupid but I had to ask.

I'm building a database of US states, cities, streets AND their basic relations to each other. I already have state and city tables.

Where can I find street data to purchase?

I know nothing about GIS and I'm too busy to learn it. So I have to buy it i believe.


texts: Editing shapefile adding fields with values with pyshp

texts:
I'm trying to add some fields to a shapefile using pyshp. The code I have to far is:
import shapefile

# Read in our existing shapefile
r = shapefile.Reader(r"C:\...\Edit_SHP\TracksLines.shp")

# Create a new shapefile
w = shapefile.Writer(r"C:\...\Edit_SHP\TracksLines_edited.shp")

# Copy over the existing fields
w.fields = list(r.fields)

# Copy over the records
for shaperec in r.iterShapeRecords():
 w.record(*shaperec.record)
 w.shape(shaperec.shape)

w.field('test')

The issue is that when I try to add a new field all the records already present on the shape are messed up.
How can I keep the same numbers of records, when adding new fields?


texts: Always getting "unable to geocode" in Nominatim?

texts:
I've installed Nominatim and everything seems to work fine, but whenever I try to use the reverse geocoding API, all I get is "unable to geocode". Strangely, the search function works well and returns lots of information for a location.

Example: The search for the coordinates 40.724362, -74.002339 correctly returns an address in NYC when entered in the search field:

But when I run

http://localhost/reverse.php?format=json&lat=40.724362&lon=-74.002339&zoom=18&addressdetails=1

The response is:

{"error":"Unable to geocode"}

What is the problem?


texts: What are PyQGIS Functions for MakeFeatureLayer and SelectLayerByLocation?

texts:
I have a script that I had created using ModelBuilder, but our field personnel use QGIS 2.0.  Below is a snippet of the code. I am wanting to know if there is an equivalent using pyqgis so that I can add it to my plugin.

# Process: Make Asset Feature Layer
arcpy.MakeFeatureLayer_management(asset_wps, Asset_Selection, "\"ASSET_TYPE\" = 'Secondary Ped' OR \"ASSET_TYPE\" = 'XFMR'")

# Process: Make PC Feature Layer
arcpy.MakeFeatureLayer_management(Proposed_Conductor, PC_Selection, "\"CONDUCTOR_\" = 'Underground Secondary' OR \"CONDUCTOR_\" = 'Underground Service'")

# Process: Select Layer By Location
arcpy.SelectLayerByLocation_management(Asset_Selection, "INTERSECT", PC_Selection, "1.5", "NEW_SELECTION")

matchcount = int(arcpy.GetCount_management(Output).getOutput(0))
print('E1012: {0}').format(matchcount)

So what I am doing here is just making a selection from two features and then checking where they intersect with that selection and finally getting a count.

How can I replace MakeFeatureLayer and SelectLayerByLocation using PyQGIS?


texts: What format is lat/long stored in OSM PostGIS?

texts:
Question:
How do I to get the lat, long for all the toll booths stored in OSM for Virginia?

Current Query:

> gis=# select * from planet_osm_nodes where ('toll_booth' = ANY(tags));

     id     |    lat    |    lon     |                            tags
------------+-----------+------------+------------------------------------------------------------
  496668427 | 448545276 | -881565606 | {name,"Airport Parking",barrier,toll_booth}
  525372041 | 463020204 | -874380697 | {name,"Swift Run Gap Entrance Station",barrier,toll_booth}
  526241181 | 470778112 | -870435097 | {name,"Front Royal Entrance Station",barrier,toll_booth}
  534508045 | 458513289 | -877777774 | {name,"Rockfish Gap Entrance Station",barrier,toll_booth}
  563894077 | 470183057 | -857852817 | {barrier,toll_booth}
  ...

Rest of the table

Additional Info

OSM Import Log
Failed Attempts
Tutorial -- manually building a tile server
Data source -- Cloudmade/virginia.osm.bz2.


texts: How to visualize the vertices of a graph in QGIS?

texts:
I tried to visualize the graph generated by osm2po and rendered that graph in qgis. Now what I am trying to do is plot the vertices in the graph. So that I can visually see where the nodes/vertices are. I have attached a screenshot of the visualization I have up to now.

I want to actually see the nodes in this

I followed the suggestion as given by underdark. However, I could see that the number of nodes and edges is less than that in the database. In the database the number of nodes = 32 and the number of edges = 39. However in my case it is 31 nodes only and 38 edges only. What could be the reason?


texts:
I am trying to align two raster grids in R. Once aligned I would like to be able to add them together.

I have tried to check whether making a stack would work:

grid_snap <- stack(habi_sdw, Pop_sdw)

And I get the following error:

  Error in compareRaster(x) : different extent

The raster grids have the following properties:

show(habi_sdw)
# class       : RasterLayer
# dimensions  : 9187, 9717, 89270079  (nrow, ncol, ncell)
# resolution  : 0.00892857, 0.00892857  (x, y)
# extent      : -28.83706, 57.92186, -36.02464, 46.00214  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
# data source : C:\Users\di39\AppData\Local\Temp\R_raster_di39\raster_tmp_2015-08-12_172902_12860_17067.grd
# names       : layer
# values      : 0, 333707.6  (min, max)

show(Pop_sdw)
# class       : RasterLayer
# dimensions  : 10143, 8858, 89846694  (nrow, ncol, ncell)
# resolution  : 0.008333333, 0.008333333  (x, y)
# extent      : -17.53524, 56.28143, -46.97893, 37.54607  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
# data source : C:\Users\di39\AppData\Local\Temp\R_raster_di39\raster_tmp_2015-08-12_170421_12860_12760.grd
# names       : pop2010ppp
# values      : 0, 128925.9  (min, max)

Using alignExtent() in the raster package seems not to be the correct approach.

Do I need to resample because the resolutions are slightly different?

(0.00892857 x 0.00892857) vs (0.008333333 vs 0.008333333)


texts: Visualize a polyline with decreasing opacity towards its ends in QGIS

texts:
For a mapping project I am looking for a method how to visualize a polyline where the ends increase in opacity such that it will look like its fades away using QGIS (3.10.10). See the red circle in the picture below.
In essence I am looking to style a polyline the way you would a DEM; with a color ramp that has 0% opacity at the lowest and highest values but instead of measures of height it would be of length. So, like a shape burst fill but for lines, I can imagine I have to play around with the geometry generator.
One solution would be to divide the polyline in several segments at fixed distances and give each segment an increased transparency, however I am looking for a less cumbersome and more robust method.
The result is purely for visual output e.g. the attribute table is not of importance.


texts: Planet SkySat - Differences between image strips

texts:
I'm trying to understand why i'm seeing variations in the analytic (surface reflectance) imagery brightness values in Planet SkySat imagery between image strips that were collected around the same time.
I'm aware of the differences between analytic and visual products (i.e. this question), but i'm hoping someone can shed some light on why the analytic product contains these differences. Surely if it is being used for scientific analysis, the brightness values should be comparable between image strips? Example below (i've mosaicked these images into one):

Am I missing a processing step here? The visual product seems to show similar differences. My use case is image classification, so ideally I need these mosaics to contain comparable brightness values across image strips to allow me to apply a classification model consistently across the whole area.


texts:
I have just transferred my Postgis database onto an RDS server, and there is a bit of a performance hit when using the Identify tool in QGIS (3) compared to localhost. I presume this is mostly network latency, however I wonder if it could be not picking up the spatial index or something. There is a GIST index set up on the table, but identify queries are still taking around 10 seconds, which seems too long, considering rendering is considerably quicker than that.

I did look at the code here and it does seem on a quick glance that things are being done client side (and presumably calling back to the server for the data hence the delay) where it would be quicker to call the server and ask it to identify the object using the indexes.

Are there any options I can tweak that anyone knows of or is this one for the QGIS developers?


texts:
I have more experience in "traditional" (Java EE) webapp development, and I am trying to figure out a decent design for a mapping web application.

My application has a map, and users can see and edit some features on the map. Let's say they can add and modify simple rectangular areas on the map, and edit some data fields for each area. Map will most likely be implemented using OpenLayers.

I will have a relational database that will store non-map-related data. Let's say the area on the map is a "school district" with direct fields such as "number of students". Then there's another page where you can input information about teachers for that area.
How do these options compare?:

A. WFS-T is not used. All the data is stored in the traditional relational DB. When user uses the map, the base map images are read from the map server but the data about the user editable features is read from the DB. I read the coordinates for each area from my backend and draw the areas using OpenLayers tools for drawing e.g. vectors on the map. All my modifiable data (map-related or otherwise) is in the same relational DB.

TABLE SCHOOL_DISTRICT
(ID
AREA_X1
AREA_Y1
AREA_X2
AREA_Y2
NAME
NUMBER_OF_STUDENTS);

TABLE TEACHER(
ID
FIRST_NAME
LAST_NAME
SCHOOL_DISTRICT_ID);

B. Use WFS-T for storing feature data. When user uses the map, both the base map images and the editable areas are read from the map server. Now it becomes a question about where to draw the line between the mapserver/WFS-T and the traditional relational db. All the data is not suitable for storing in the map server. The "school district" data may be linked to a complex graph of tables which need to be e.g. reported on. Should I store just the geometry data (coordinates of the areas) in mapserver/WFS-T and everything else in traditional db tables? Or should the immediate data about school districts be stored in mapserver/WFS-T and rest in traditional db tables? How do these two combine? Or, can mapsers/WFS-T utilise the same tables/schema as the traditional side of the application...?

GEO_TABLE_IN_GEOSERVER SCHOOL_DISTRICT_AREA
(ID
AREA_GEOMETRY);

TABLE SCHOOL_DISTRICT_AREA
(ID
AREA_GEOMETRY_ID??
NAME
NUMBER_OF_STUDENTS);

TABLE TEACHER(
ID
FIRST_NAME
LAST_NAME
SCHOOL_DISTRICT_ID);

EDIT: adding an image to hopefully illustrate what I am talking about

We have a web application that has both

spatial data viewed and edited in an OpenLayers map interface (left side, for example geometries of "school districts")
non-spatial data which is viewed and edited elsewhere in the application (right side, for example "teacher" data)

Basically, spatial data is edited using WFS-T. The map server (GeoServer) offers WFS-T interface and ultimately talks with RDBMS using SQL. Non-spatial data is edited with means from "standard" web applications: forms, http requests, application server (such as JBoss etc) handles the requests and ultimately talks with RDBMS using SQL. Map server is not involved in this.

Both spatial and non-spatial data can reside in the same RDBMS.

What I am still unsure of is:

What data should be updated using WFS-T? Is it just the pure geometry (shapes of "school districts") or also data directly associated with the geometries ("school district's number of students -property)? Or can we use WFS-T for all updates, also non-spatial data ("teacher" data - although this does not sound nice)
If the stored geometries are simple (points in the map, or rectangular areas), is it smart to involve WFS-T and it's added complexity, or would it be sensible alternative to persist those geometries using the right side of the diagram? You would just have simple X1,Y1,X2,Y2 fields in a non-spatial table "school district", and handle storing that data in your regular web application just like you store all other data.


texts: Determining which feature is identified on ArcMap document using ArcPy?

texts:
I have to find out which feature is identified on the map, thus subsequently read the attributes of the feature, and then do some logic to plot a time series graph.
Is it possible to find out which feature is identified on the map using the identify tool, or alternatively is it possible to find a graphic below a mouse click using python in ArcGIS.


texts: Create new shapefile from Geometry generator

texts:
I have a line shape and used QGIS geometry generator to modify the line. Now I would like to create a new shape file from the new style. Is that possible in QGIS?

I'm using QGIS 3.4


texts: Export region polygon from openstreet map

texts:
Could someone inform whether it is possible to export the region Loddon Mallee from openstreetmap to geojson format?

http://www.openstreetmap.org/relation/4246154#map=7/-35.791/143.064


texts: What happened to the "Enclosing Zip" field in the ArcGIS Zip Code Data?

texts:
In the "Data and Maps" data that is included with ArcGIS Desktop, the zip code points file, 'zip_usa' (or for example, the layers referenced here: http://www.arcgis.com/home/item.html?id=1eeaf4bb41314febb990e2e96f7178df ) used to have a field called "ENC_ZIP".

For point zip codes (such as P.O. Boxes or single buildings that have their own zip code), this told you what zip code the point zip could be found inside.

When I look at current versions of this data, that field has been removed.

Does anyone know where this data went?  Is there any other good source for this information?


texts:
I downloaded several Landsat images from some years between 1985 and 2014, all of them with good quality. Nevertheless the period from 2003 to 2013 is not worth it. I checked all images. Most of the images have black strips or color distortions.

Checking on the EarthExplorer page I get the same result. Is there something wrong with that period of detection, do I just have bad luck with my area, or am I doing something wrong?


texts:
Is there a function that produces an ISODATA unsupervised classification in QGIS?

The only unsupervised classification function I can find is a KMeans.


texts:
The basic issue is that I have shapefiles saved in WGS84/Pseudo Mercator (EPSG: 3857) and some DEM data saved in NAD83/UTM Zone 10N (ESPG: 26910), but can't get them to line up.

I'm using QGIS and have my file CRS set to WGS84 (3857) with on-the-fly projections enabled. My global settings allow on-the-fly to automatically be applied to any new layers added that aren't in the drawing CRS.

The DEM files that I'm having trouble with are from the National Elevation Dataset (NED) and were obtained through the NRCS Geospatial Data Gateway. I've tried to change the DEM layer CRS in the layer properties. I've also tried to change the CRS by right-clicking the layer and choosing "Set Layer CRS" and I've tried to "save as..." and alter the CRS for a new layer file. None of those have worked.

Then I tried to use GDALWARP to reproject the DEM from NAD83 to WGS86, but the DEM just ends up in exactly the same spot.

Trying to change the CRS of the drawing and the shapefiles to meet the DEM is similarly counter productive. They are outside the bounds of the NAD83 UTM Zone 10N, so they just disappear.

I'm pretty new to QGIS and GIS in general, so it's likely that I'm making a real rookie mistake here. I'm using QGIS 2.8.1 Wien on a Windows 8 computer if that helps.


texts: Remove some properties after inner Join in Google Earth Engine

texts:
I want to join two featureCollections and remove the ID field in order to run linearRegression in GEE. I found codes to clean the featureCollection after innerJoin, and snippet to remove properties from FC. However, when I try on my own data, all the properties from FC were removed!
Here is my code
// Create the primary collection.
var primaryFeatures = ee.FeatureCollection([
  ee.Feature(null, {foo: 0, label: 'a'}),
  ee.Feature(null, {foo: 1, label: 'b'}),
  ee.Feature(null, {foo: 1, label: 'c'}),
  ee.Feature(null, {foo: 2, label: 'd'}),
]);

print('primary fc',primaryFeatures);

// Create the secondary collection.
var secondaryFeatures = ee.FeatureCollection([
  ee.Feature(null, {bar: 1, label: 'e'}),
  ee.Feature(null, {bar: 1, label: 'f'}),
  ee.Feature(null, {bar: 2, label: 'g'}),
  ee.Feature(null, {bar: 3, label: 'h'}),
]);

print('2nd fc',secondaryFeatures);

// Use an equals filter to specify how the collections match.
var toyFilter = ee.Filter.equals({
  leftField: 'foo',
  rightField: 'bar'
});

// Define the join.
var innerJoin = ee.Join.inner('primary', 'secondary');

// Apply the join.
var toyJoin = innerJoin.apply(primaryFeatures, secondaryFeatures, toyFilter);

// Print the result.
print('Inner join toy example:', toyJoin);

// I tried
// clean the join
// this snippet is from stackexchange
function cleanJoin(feature){
  return ee.Feature(feature.get('primary')).copyProperties(feature.get('secondary'));
}
var Joined = toyJoin.map(cleanJoin);
print('joined',Joined)

// remove one properties
// this snippet is also from https://gis.stackexchange.com/questions/321724/removing-property-from-feature-or-featurecollection-using-google-earth-engine

var removeProperty = function(feat, property) {
  var properties = feat.propertyNames()
  print(properties);
  var selectProperties = properties.filter(ee.Filter.neq('item', property))
  return feat.select(selectProperties)
}

var removeJoined = removeProperty(Joined,'foo');
print('remove joined',removeJoined) // all properties were removed!!


texts: Adding  to Leaflet 1.0.0 layers control throws an error

texts:
My intention is to add transparency sliders to each layer in a L.control.layers. I was going to use a simple HTML5 slider <input>, demonstrated here:

http://geohackers.in/2013/10/leafelt-layer-opacity-control-using-html5-slider/

However, when I add an <input> to the legend string in a L.control.layers, I get an error:

leaflet.js:12457 Uncaught TypeError: Cannot read property 'layer' of undefined(…)
_checkDisabledLayers @ leaflet.js:12457
_addItem @ leaflet.js:12411
_update @ leaflet.js:12325
onAdd @ leaflet.js:12174
addTo @ leaflet.js:11610
(anonymous function) @ index.html:144

Removing the <input> fixes the error. What's the explanation?


texts: Static Image Layer Performance Issues - Open Layers 3

texts:
I am receiving data points in lat/lon. I need to plot all of these on a static image in OpenLayers 3. I am not using OSM or anything else, my image should be the base layer. I got it working by adding the static image layer, setting the extent to be the four corners in lat/lon and then plotting my features on top of it. However this runs in the browser very slowly, and on mobile it will crash the browser after a couple seconds.
What is the proper way to do this? Should I be using an overlay instead? Do I need to convert my lat/lon coordinates to pixels?

    var ovProj = ol.proj.get('EPSG:4326');
var myStaticImageLayer = new ol.layer.Image({
    source: new ol.source.ImageStatic({
        url: '../_images/mallSitePlan.png',
        imageExtent: [-121.90320739419553, 37.409945396290674, -121.89234011506653, 37.41962634032544],
        projection: ovProj
    })
});

var view = new ol.View({
    center: ol.proj.transform([-121.90320739419553, 37.409945396290674], 'EPSG:4326', 'EPSG:3857'),
    zoom: 18,
    enableRotation: false
});

var map = new ol.Map({
    target: 'map',
    layers: [myStaticImageLayer],
    view: view
})


texts: How to calculate the azimuth of an offset point/line?

texts:
After receiving and answer to this question, I am stuck with another problem.
For context, my end goal is to assess windward/lee side of a forested edges (line segments) in respect to a dominant wind direction. Following the answer in the linked question, I have created an offset point on the inside edge of each line segment. I now am trying to evaluate these line segments based on the azimuth of the offset (segment midpoint to offset midpoint).
How can I reference the offset point in an expression in a field calculation on the original line segment?

How do I reference offset_point if it is in another layer?
If it is another feature in the same layer?
psuedocode:
offset_azimuth = degrees(azimuth(segment_midpoint, offset_point))
if(offset_azimuth = (windDir +/- 45), "Lee", "Windward")


texts:
I am using ArcGIS desktop 10 and I would like to change a 2D polygon into 3D building.

Does anyone have any idea from where I can download the 3D building generator tool for ArcGIS Desktop?

I read about this tool here.

This tool gives various features and can build complex roof tops for building sites. With this tool I can develop a rounded roof in 3D ..etc

Moreover, any suggestions how I can import this 3D buildings in a web service using asp.net and C#?

Any help please?


texts: Saving dataframe as dbf with simpledbf in Python?

texts:
I am doing a number of joins outside of ArcGIS as pandas seems to outperform it. I read in the dbf file with simpledbf, convert it to a dataframe, then do the join:

import pandas as pd
from simpledbf import Dbf5

#Source of the join
data_source=pd.read_excel('C:\\Users\\myfolder\\mysource.xlsx')

#Target of the join
dbf = Dbf5('C:\\Users\\myfolder\\mydbf.dbf')
d_b_f = dbf.to_dataframe()

#Join
d_b_f_2=pd.merge(d_b_f, data_source, on='COLUMN', how='left')

How can I save this dataframe as a dbf and overwrite the existing dbf?

I expected something along these lines to work:

final_dbf = Dbf5('C:\\Users\\myfolder\\finaldbf.dbf')
final_dbf = d_b_f_2.to_dbf()

but an error pops up saying that 'C:\Users\myfolder\finaldbf.dbf' does not exist. What do you suggest as a workaround? I wouldn't want to look into other packages, as I am sure there must be a way of doing this within simpledbf itself.


texts:
I am trying to add multiple tables using tabs to one of my viewers but can't get it to work.  Below is the code that I have added to my viewer.js widget option (it adds a table for 'Site 2' but not 'Site 1').

attributesTable: {
            include: true,
            id: 'attributesContainer',
            type: 'domNode',
            srcNodeRef: 'attributesContainer',
            path: 'gis/dijit/AttributesTable',
            options: {
                map: true,
                mapClickMode: true,

                // use a tab container for multiple tables or
                // show only a single table
                useTabs: true,

                // used to open the sidebar after a query has completed
                sidebarID: 'sidebarBottom',

                // optional tables to load when the widget is first instantiated
                tables: [
                    {
                        title: 'Site 1',
                        topicID: 'site1Query',
                        queryOptions: {
                            queryParameters: {
                                url: 'https://your.rest.url.com/gis/rest/services/XXXX/YourMapService/MapServer/9',
                                maxAllowableOffset: 100
                                },
                            idProperty: 'OBJECTID'

                        },

                        title: 'Site 2',
                        topicID: 'site2Query',
                        queryOptions: {
                            queryParameters: {
                                url: 'https://your.rest.url.com/gis/rest/services/XXXX/YourMapService/MapServer/10',
                                maxAllowableOffset: 100
                               },
                            idProperty: 'OBJECTID'

                        }
                    }
                ]
            }
        },


texts: how can I obtain the areas where my rasters overlap, using Qgis?

texts: I have a bunch of raster files showing the distribution of different species of mammals, and I would like to obtain the areas where they overlap, the areas in whitch more than one species can be found


texts:
From an ArcGIS Server (9.3.x or 10.x) instance, can anyone confirm if REST HTTP response headers are modified when the data within a layer is updated?

For example, if I query a layer at time x and get a 200 OK response and E-Tag y and I query the same dataset at a later date, I'll generally receive a 304 Not Modified response with the same E-Tag. However I'd expect the E-Tag to change and to receive a 200 OK status if the layer data changed, is this the case?

If this is not the case, is there a way to test if data has changed since a particular date without actually pulling the data down.


texts: Opening Esri Personal Geodatabase (*.mdb) using QGIS?

texts:
Is it possible to open an Esri Personal Geodatabase (*.mdb; ArcGIS 10.2) in QGIS 2.6?
Previous versions were right.


texts: Converting branched multiline into series of single lines in shapely

texts:
If I have a multiline string representing a branched network with a common starting point (red network in pic below):

Does anyone know how to convert this to a set of individual single linestrings in shapely? With each individual linestring starting at the common starting point and ending at each specific end point? Assuming I have a shapely multiline object as a starting point.


texts: Random Stratified Sample in QGIS and Scale

texts:
I need to sample 7 raster layers all 30 m resolution (as some of the remotely sensed products I'm using are Landsat).

I have all of my layers ready to go, but I am not sure how I can conduct a stratified random sample. The forest I'm interested in has varying forest ages (primarily old) with fewer amounts of younger forests that as shapefiles look like hundreds of small polygons.

Additionally, I have 3 types of geology for the region. If I generate random points using the Random Points tool, I've noticed that I end up with lots of points and may be sampling a pixel more than once, which seems incorrect.

Is there a way for me to sample my raster files. Shouldn't my points be the same size as my pixels? These data will be put into a model I'm building, so having representative samples for these different classes is important.


texts: Best way to store, elaborate and query raster data

texts:
I have created multiple raster images with GRASS-GIS (~150Gb), that I'll need to multiply for a certain coefficient, and query the results that need to by saved just temporally. The portions of the raster image that need to be queried are identified by the selection of some polygons contained in a multiple polygon vector. The queried data then need to be elaborated by a simple python script for each of the selected polygons and return the result in a geojson format.

Which is the best way and tools for this operation?


texts: Displaying animated data over time with the Google Maps API and JS

texts:
I'm building an app using Google Maps API and AngularJS.
I've managed to display data both as pointers, clusters and heat maps and to reload the maps in real time when the data changes.

Now I was thinking about the possibility to actually display a dynamic/animated heatmap that, for example, given a certain array of objects with with date property, displays a heatmap that changes over time, so you see a progression.

After many google searches I haven't found anything built with Google Maps that does something similar.

Could anyone tell me if it is even possible? And if not, which technology would you recommend for such a map?


texts:
I wonder is it possible to get lowercase GUID in attribute table for my feature class and how to do that?
So, when I calculate appropriate field (type - guid) with a following script:
def ID():
  import uuid
  return '{' + str(uuid.uuid4()) + '}'

I get uppercase string, like that: {CC82E261-1B32-4189-A22F-4FD33C0A5872}. But I need {cc82e261-1b32-4189-a22f-4fd33c0a5872}.
I tried to apply function !globalid!.lower() (works only in editing mode), but when I try to save it, ArcMap returns uppercase string.


texts: How to show two layers with different projections in QGIS

texts:
I have a shapefile that define X and Y in meters, not lat and lon. Then I have another file (CSV) that stores points locations as lat and lon, e.g. 30.4, 2.2.

I want to load these two files as two layers in QGIS 2.10. First I added a shapefile as a Vector Layer. Then I added a CSV file as a Delimited Text Layer.

I can only see a map (shapefile). How can I correctly project lat and lon? I tried to set different CRS, e.g. EPSG: 3857, for the Delimited Text Layer, but nothing works.


texts: Converting DWG surface to line feature class in file geodatabase using FME?

texts:
I have a DWG file with autocad_surfaces.

I want to export these to a line, so that I can use it in ArcGIS Pro.

How can I convert the features from surfaces to lines in an ArcGIS File Geodatabase using FME 2018?


texts: How to open MDB personal geodatabase in JAVA

texts:
I'm trying to open a MDB personal geodatabase using the ArcGISSDK java10.2.4 doing the following:

String pathMDB="C:\\Documentos\\NetBeansProjects\\Pruebas\\OpenPersonalGeodatabaseMDB\\GDB_CI6_L01.mdb";
        Geodatabase mdbGeodatabase=null;
        try {
                mdbGeodatabase= new Geodatabase(pathMDB);
        } catch (FileNotFoundException ex) {
            Logger.getLogger(OpenPersonalGeodatabaseMDB.class.getName()).log(Level.SEVERE, null, ex);
        }

I get the error "File opened that is not a database file". I tried with two diferents files whith the same result. Is there another way to do it ?.


texts: Identifying lines of one feature class that touch another feature class with their endpoints in ArcMap?

texts:
I want to find out and select the features of one feature class, which both endpoints touch another feature class.

In the picture above it would be the left pink line which touches feature of another feature class on both ends.

The best would be, if there is a way to find out the right features automatically.
This question ( Distinguishing lines that intersect from lines that touch? ) could also be a solution, but I have difficulties using the Python Script. I want to use this as a script tool but I don't know how to select parameters. I thought they will be added automatically after adding the Python script, but it didn't work out.

How can I do this?


texts:
I want to use Finland national WMS map (MML Kapsi) as background map in FME inspector. Is it possible? What information is needed?

WMS URL is: http://tiles.kartat.kapsi.fi/peruskartta


texts: Using intersect method of ItopologicalOperator gives error 0x80040239?

texts:
I would like to intersect a polygon and a polyline to know if the polyline is in the polygon.

The code works but when i do it on more than 100 000 times in a for, there is this error :

System.Runtime.InteropServices.COMException (0x80040239): Exception de HRESULT : 0x80040239
    à ESRI.ArcGIS.Geometry.ITopologicalOperator.Intersect(IGeometry other, esriGeometryDimension resultDimension)

My code :

ITopologicalOperator4 topoOperateur = secteurPsr.Shape as ITopologicalOperator4;
topoOperateur.IsKnownSimple_2 = false;
topoOperateur.Simplify();
IGeometry geomPsr = currentPsr.Value.Shape;
IPolyline polylinePsr = geomPsr as IPolyline;
(geomPsr as ITopologicalOperator2).Simplify();
IGeometry geomResult = null;
try
{
    geomResult = topoOperateur.Intersect(geomPsr,esriGeometryDimension.esriGeometry1Dimension);
 }


texts:
I have created my own repository for QGIS. If I now install my plugin from this repository, I get the following error message:

  Plugin has disappeared: The plugin seems to have been installed but I don't know where. Probably the plugin package contained a wrong named directory. Please search the list of installed plugins. I'm nearly sure you'll find the plugin there, but I just can't determine which of them it is. It also means that I won't be able to determine if this plugin is installed and inform you about available updates. However the plugin may work. Please contact the plugin author and submit this issue.

When I restart QGIS and activate the plugin, it works without problems.

My repository looks like this:

<?xml version="1.0"?>
<plugins>
    <pyqgis_plugin name="Name" version="0.1">
        <description>xxx</description>
        <homepage>www.xxx.com</homepage>
        <qgis_minimum_version>3.0</qgis_minimum_version>
        <file_name>plugin.zip</file_name>
        <author_name>Name</author_name>
        <download_url>http://xxx.com/plugin.zip</download_url>
    </pyqgis_plugin>
</plugins>


texts: gdalbuildvrt and gdal_translate. Incorrect positioning of the tif

texts:
My program builds layers from TAB and TIF files, using gdalbuildvrt and gdal_translate tool. The typical TAB file looks like so:

!table
!version 300
!charset WindowsCyrillic

Definition Table
  File "12-a.tif"
  Type "RASTER"
  (-2,14) (186,180) Label "Point 1",
  (-2,13) (211,6072) Label "Point 2",
  (-1,13) (6108,6052) Label "Point 3",
  (-1,14) (6088,161) Label "Point 4"
  CoordSys NonEarth Units "km"
  Units "km"
  RasterStyle 1 86
  RasterStyle 2 100
  RasterStyle 4 1
  RasterStyle 5 0

gdalbuildvrt creates mosaic.vrt file and gdal_translate creates a georeferenced tif. This is how the commands look like:

$ gdalbuildvrt -input_file_list files.txt mosaic.vrt
$ gdal_translate -co "TILED=YES" -co "BIGTIFF=YES" -co "NBITS=1" -co "COMPRESS=CCITTFAX4" mosaic.vrt layer.tiff

This is for example a part of mosaic.vrt file, which corresponds to the above 12-a.tif:

<SimpleSource>
  <SourceFilename relativeToVRT="0">12-A.tif</SourceFilename>
  <SourceBand>1</SourceBand>
  <SourceProperties RasterXSize="5990" RasterYSize="5995" DataType="Byte" BlockXSize="5990" BlockYSize="1" />
  <SrcRect xOff="0" yOff="0" xSize="5990" ySize="5995" />
  <DstRect xOff="22.3535484203381" yOff="0" xSize="5997.59263406113" ySize="6002.59897181911" />
</SimpleSource>

The problem is, when I publish my layer, I see that the tif (for example, this very 12-A.tif), is incorrectly positioned (about 20 meters offset along y axis and about 70 meters offset along x axis). One important thing I should add, is that my algorithm does some conversion of the initial TAB files - first of all, it converts initial local coordinates (with unknown SRID) to some well know coordinate system and it converts absolutely correctly (I checked the precision of coordinates dozens of times and they are perfect) and it switches from "km" units to "m" units. All the rest inside TAB files is left untouched, including pixel offsets like in this line "(-2,14) (186,180) Label "Point 1"". Here is how it looks like:

The gray point on the picture marks where the top left corner should be and, as you can see, it is far away - about 20 meters to the north and 70 meters to the east.

So, my question is what may be wrong with that and how can I possibly fix it?

EDIT

If I apply gdalinfo to my 12-a.tif, I see that coordinates of Upper Left corner are different from those in TAB file. So, is it a correct behaviour of GDAL?


texts: How to workaround truncation of text description for elements in qgis-modeler 2.18.15?

texts:
I am trying to make a QGIS processing workflow through the Graphical Modeler (QGIS 2.18.15). I have completed it, and wish to export it as a *.PNG. However, the elements and algorithms I added to the modeller are visually incomplete (unable to see the whole text written per element).

At present, it looks like this

Since I wish to export it, I want it to visually look like the following example sourced online.

How do I customize it such that my text per element is not truncated?
And, is it possible to customize the aesthetics of the QGIS graphical modeller (within our outside of QGIS)?


texts: ArcGIS Server SOC Cannot Connect To Data

texts:
I've been successfully publishing services to AGS.  SOC is set up to view all the data that is used in published services.  All of a sudden, I can no longer publish a service and have recieved the following error.

  Possible cause of connection failure: the soc cannot resolve the path to the data, or the soc account has not been granted read permissions of the data. In the case of ArcSDE, direct connect may not be configured properly, or the soc account may not have permission to access the database.

These are some things I have noticed:

If I drag and drop data into a new map document from an existing successfully published map, it works.  If I add that same dataset, but from Catalog, it does not work.
All the data from a successfully published map is referenced via a UNC Path and show up as medium warnings when analyzed through the publishing tool.  If I add a dataset from Catalog (or Add Data tool) and analyze the map, it indicates that it is not referenced via a UNC path.
If I create a map document that contains a dataset (again, dragged and dropped from a successuflly published map) that is UNC referenced and a dataset add through Catalog/Add Data Tool, it will publish but will leave out the dataset that is not UNCed.

From the image below, the three layers come from the same GDB.  The top layer was dragged and dropped from another map document and will show up in the published service.  The bottom layer is the same dataset except it is add to the map through the Add Data tool.  This is where the above error is coming from.

All data that is brought in from Catalog/Add Datat tool is no longer UNC referenced.

Looks like something broke, any thoughts? Maybe reboot AGS?


texts:
How would you implement geographic profiling in QGIS?

There doesn't seem to be any pre-existing plugins and my Python skills are not quite good enough to implement it myself. Specifically, Rossmo's Formula would be a good contender for creating an accurate crime profile map.


texts: How not to lose columns in a result of a pixel extraction, if some column is equal to zero?

texts:
I have reclassified some classes of a raster into two classes (natural and anthropic) and I have to extract the pixel proportion of each class for some study areas that I have. However, there is a time when, among all my study areas, one of them did not have values for the anthropic class (it's ok, this is normal), but, the code I am using excludes this column and I don' t want to exclude it, because I will bind (rbind) with results from my others study areas, and they should have the same length ( columns number).

#The data raster
map<-raster("prodes.tif")

raster structure

study_area<- readOGR("study_area)

Shapefile structure

 #extract the pixels values to my study area
 ( extrc <- extract(map, study_area) )

      #here calculted the prportion of each class (anthropic and natural)
      ( class.prop <- lapply(extrc,FUN = function(x) { prop.table(table(x))}))

  # now I creted use this function to tranform  "class.prop" in a Data Frame
   rbind.fil <- function(x) {
   nam <- sapply(x, names)
   unam <- unique(unlist(nam))
   len <- sapply(x, length)
   out <- vector("list", length(len))
   for (s in seq_along(len)) {
   out[[s]] <- unname(x[[s]])[match(unam, nam[[s]])]
    }
   setNames(as.data.frame(do.call(rbind, out), stringsAsFactors=FALSE), unam)
  }

  ( p.prop <- rbind.fil(class.prop) )

      # I inser a new IDcolumn
      p.prop$ID<-seq_along(p.prop[,1])

      #change the classes names
      names(p.prop)[names(p.prop) == '0'] <- 'anthropic'
      names(p.prop)[names(p.prop) == '1'] <- 'Natural'

It's missing the anthropic column

      #Do the same to my study area shape
      study-area$ID<- 1:length(study-area)

      freq <- merge(study-area, p.prop) #add to polygons

      frequenc<-as.data.frame(freq) # transfor to data frame

      results <- rbind(results, frequenc) #bind with my previous results

      write.table(results, file="resultados.txt", sep="\t", row.names=F) #save as txt


texts: How to identify the direction of travel for ways without "oneway" tag

texts:
I am using Overpass API to download ways and nodes information for my application and I have realized that not all the ways have information related to the direction of travel.

How can I identify if the way supports direction in both sides/with line direction or against direction when no tags are attached?

Example url: https://lz4.overpass-api.de/api/interpreter?data=[out:json];way(50.05288170227778,8.69143545627594,50.05597466228186,8.694798946380615);out%20meta;


texts:
I have computed some isochrones using OSM route plugin but when I am trying to use the plugin again it gives me an error: An error has occured while executing Python code.

Is it working for you?


texts: Rendering geoJson polygons in google maps android apiv2

texts:
I'm in the process of publishing spatial data stored in a postgis database to a google maps based android app and am trying to put together a plan for rendering my data in the most efficient manner possible.

The portion of my dataset I'm concerned about consists of about 200,000 polygons that will be filtered by the user name, bounding box etc.  I can see situations where I would like to display up to 1000 heavily simplified polygons when zoomed out, but I doubt that I will be displaying more than 20 on a screen without simplifying the data.

My current strategy is to query my database based on the filtering factors mentioned above, return the query as a geoJson and render the resulting data directly into google maps.

Is this an appropriate strategy or should I be looking at using another service to help render my data before it is loaded into google maps.  If it is an appropriate strategy I could use some advice on how to render geoJson polygons in the google maps android api.

I'm quite comfortable with the SQL end of this project but I'm not as experienced with javascript or the google maps android apiv2.

Thank you for the help.


texts: How to plot a list of points on Overpass turbo (OpenStreetMap)

texts:
I am using overpass turbo web http://overpass-turbo.eu/#

after typing in

[out:csv(::lat,::lon;false)];
relation(2081626);>;out;

I get already a list of coordinates, under "data" tab

48.0786156  11.5510212
48.0769149  11.5502003
48.0763526  11.5505930
48.0768127  11.5502292
48.0761811  11.5499233

...

How could I plot this list of coordinates on "Map" tab as a route, on the map


texts: add vb6 tool to arcmap 9.2 (on windows 7 64bit machine)

texts:
I have a vb6 project which develop a tool for use in arcmap several years ago. And I managed to add the tool in arcmap 9.2 (windows xp) by:

In arcmap -> Customize -> Add from file -> choose the dll file -> tool be added successfully.

But when I follow the same steps to add the tool in arcmap 9.2 (windows 7 64bit), I get "No new objects added".

I want to check is there a way to add a vb6 tool to arcmap 9.2 on windows 7 64bit machine?


texts: Present change in cumulative sum of point in polygon by month

texts:
I have a data table that contains points in a city, and a polygon for statistical areas within the city. For each point, I have the date it first appeared, from 2009 to 2016.
I want to create in some way a map that presents the cumulative sum of points for each statistical area, by month. Ideally, it would show the growth of number of points for each area with color change.

edit:
I use QGIS 2.18. The database is exported from R to csv, and then imported to QGIS. The time format is dd/mm/yyyy though it won't be difficult to reformat it if necessary.


texts: Controling polygons centroids alignment in QGIS?

texts:
I am using Polygons Centroids feature in QGIS and it does an amazing job. My goal though is to set the points on the left side of the polygon and not in the center. Is there a way to control the alignment of this feature?
The polygons look like those in the images. The polygons come in different angles.

In the following, I would line the points to be on the left:


texts: Retrieving ArcGIS Collector points when service definition was deleted?

texts:
I made a service to collect data points in the field, but then right before leaving updated the fields and deleted the old service, and published my new service. I assumed the Collector for ArcGIS app would automatically update. I am now stuck with data that belongs to a service definition that no longer exists.

Is there any way that I can upload these points without the original service definition they came from?


texts: GeoServer - WMS v WFS (polygon label issue)

texts:
I am hoping that when I switch over to WFS that my multi-polygon labelling issue will disappear.

Currently I am using WMS, and I believe the reason I am seeing multiple labels within my polygons is because each tile is being rendered in its own silo. Therefore the server thinks that there needs to be a label for each tile inside a given polygon.

I know there are the centroid solutions, except that solution does not work for me since I want the labels inside the polygon and not in space outside a non-symmetric polygon.

Will WFS solve this issue?


texts: ArcGIS JS - Missing "style" when printing dynamic layer

texts:
Im having an issue trying to get the url of the map being printed with a dynamic layer. WIth just tiled services and basemap it prints just fine.

 // USED FOR PRINTING OF MAPS
    function print(map_div_id)
    {
        var printUrl = "http://prodgis.agriculture.purdue.edu:6080/arcgis/rest/services/Utilities/PrintingTools/GPServer/Export%20Web%20Map%20Task/execute";
    printTask = new esri.tasks.PrintTask(printUrl);

    var template = new esri.tasks.PrintTemplate();

    template.format = "JPG";
    template.layout = "MAP_ONLY";
    template.preserveScale = true;

    var params = new esri.tasks.PrintParameters();
    params.map = map[map_div_id];
    params.template = template;

    printTask.execute(params, printComplete);
}

My Error I get is

0: "Error executing tool.: Layer "prime_soil_map_graphics": Missing 'style' in symbol.↵Failed to execute (Export Web Map).↵Failed to execute (Export Web Map Task)."

Any ideas what is happening ?

Here is a JS fiddle to test with just press the print button !

http://jsfiddle.net/YMzfM/1/


texts:
I have created the following function that measures the distance from one point in a GDF to all points in another GDF and returns a table back with the shortest distance for each point. It works well for one point however I neglected the fact I have a table of 4000 points and so it takes 10 mins. I have ran it in PostGIS and can get it down to less than a second. Is there a way to do this in Python that could match the PostGIS speed?
def get_distance_to(gdf_in, aoi_df, aoi):
   dist_df_list = list()
   for row in range(len(gdf_in)):
      single_row = gdf_in.iloc[row]
      distances = aoi_df.geometry.distance(single_row.geometry)
      dist_list = distances.to_list()
      closest_aoi = min(dist_list)
      single_row["dist_to_"+aoi] = closest_aoi
      df = single_row.to_frame().T
      dist_df_list.append(df)
completed_distances = pd.concat(dist_df_list, ignore_index=True, sort=False)
return completed_distances

my input tables looks something like this

and the output table looks like this


texts: leaflet-kml loads entire KML of 15-20 MB at client browser in one go or in chunks as and when user drags to other locations of country?

texts:
I am trying Leaflet KML plugin as below:
https://github.com/windycom/leaflet-kml
Here is full code:

<meta name="viewport" content="width=device-width, initial-scale=1.0,  maximum-scale=1.0, user-scalable=no" />
    <link rel="stylesheet" href="http://unpkg.com/leaflet@1.4.0/dist/leaflet.css" />

I have KML file of my whole country that size around 15-20 MB. As the code shows the way I loaded KML, I have a doubt that the map gets downloaded and load full KML of 15-20 MB on user's browser side in one go or the KML loads in chunks as and when user drags to other locations of country. I don't want KML to be loaded entirely in one go and slow down loading at browser side. Is it possible to load in chunks as and when user drags map to other locations of country?


texts: Convert several SpatialPointsDataFrames to raster

texts:
There are 100 SpatialPointsDataFrames in my Workspace. I loaded them on this way:

filenames <- list.files(path="",
                        pattern="XYhectareTravelTimes_ez+.*shp")

for(i in filenames){
              filepath <- file.path("/",i)
              assign(i, readShapePoints(filepath))

They are called like this:

XYhectareTravelTimes_ez10.*shp
XYhectareTravelTimes_ez11.*shp
XYhectareTravelTimes_ez12.*shp

etc.

How can I convert them to rasters looping through the Workspace?

I am very new R user and hoping to find some help.
thanks a lot.


texts: Adding layers using ArcGISDynamicMapServiceLayer

texts:
I am right now adding individual layers using FeatureLayer but it takes quite some time to load as the map is for US Counties (3076 count). I tried to load it using ArcGISDynamicMapServiceLayer but I was not successful. It seems that ArcGISDynamicMapServiceLayer can load mapserver but not the layers in the mapserver.

I used the following code but it did not work.

var ftrLyr = new esri.layers.ArcGISDynamicMapServiceLayer(url, {
    visibleLayers: [1]
});

Is it not possible to load layers in a mapserver using DynamicMapService?


texts: How to output geoJSON file from mysql database shapefile?

texts:
I have a mysql spatial database ('mydb') of one table ('tst'), I mean by spatial. I need to output the geoJSON file from this spatial mysql table. I downloaded geoPHP library and also MySQL to GeoJSON script which I adapted to my database settings, like below:
<?php
/**
* Title:   MySQL to GeoJSON (Requires https://github.com/phayes/geoPHP)
* Notes:   Query a MySQL table or view and return the results in GeoJSON format, suitable for use in OpenLayers, Leaflet, etc.
* Author:  Bryan R. McBride, GISP
* Contact: bryanmcbride.com
* GitHub:  https://github.com/bmcbride/PHP-Database-GeoJSON
*/

# Include required geoPHP library and define wkb_to_json function
include_once('geoPHP/geoPHP.inc');
function wkb_to_json($wkb) {
    $geom = geoPHP::load($wkb,'wkb');
    return $geom->out('json');
}

# Connect to MySQL database
$conn = new PDO('mysql:host=localhost;dbname=mydb','root','test123');

# Build SQL SELECT statement and return the geometry as a WKB element
$sql = 'SELECT *, AsWKB(SHAPE) AS wkb FROM mytable';

# Try query or error
$rs = $conn->query($sql);
if (!$rs) {
    echo 'An SQL error occured.\n';
    exit;
}

# Build GeoJSON feature collection array
$geojson = array(
   'type'      => 'FeatureCollection',
   'features'  => array()
);

# Loop through rows to build feature arrays
while ($row = $rs->fetch(PDO::FETCH_ASSOC)) {
    $properties = $row;
    # Remove wkb and geometry fields from properties
    unset($properties['wkb']);
    unset($properties['SHAPE']);
    $feature = array(
         'type' => 'Feature',
         'geometry' => json_decode(wkb_to_json($row['wkb'])),
         'properties' => $properties
    );
    # Add feature arrays to feature collection array
    array_push($geojson['features'], $feature);
}
header('Content-type: application/json');
//write to json file
$fp = fopen('mytable_geojson.json', 'w');
fwrite($fp, json_encode($geojson));
fclose($fp);
$conn = NULL;
?>

But when I execute the code on the browser they show that:
    out('json'); }
# Connect to MySQL database
$conn = new PDO('mysql:host=localhost;dbname=mydb','root','test123');
# Build SQL SELECT statement and return the geometry as a WKB element
$sql = 'SELECT *, AsWKB(SHAPE) AS wkb FROM mytable';
# Try query or error
$rs = $conn->query($sql);
if (!$rs) { echo 'An SQL error occured.\n'; exit; }
# Build GeoJSON feature collection array
$geojson = array( 'type' => 'FeatureCollection', 'features' => array() );
# Loop through rows to build feature arrays
while ($row = $rs->fetch(PDO::FETCH_ASSOC)) { $properties = $row;
# Remove wkb and geometry fields from properties
unset($properties['wkb']);
unset($properties['SHAPE']);
 $feature = array( 'type' => 'Feature', 'geometry' => json_decode(wkb_to_json($row['wkb'])), 'properties' => $properties );
# Add feature arrays to feature collection array
array_push($geojson['features'], $feature); } header('Content-type: application/json');
//write to json file
$fp = fopen('mytable_geojson.json', 'w');
fwrite($fp, json_encode($geojson));
fclose($fp);
$conn = NULL; ?>


texts: Why is editting in QGIS producing invalid geometry in PostgreSQL database?

texts:
I am editing a boundary in QGIS(3.10.8-A Coruña) and facing an issue when writing the back to PostGIS (10.13).
The workflow I am following is:

Download data from PostGIS connection in QGIS
Edit a boundary by: toggle editor on, switch on vertex tool and switch on snapping tool
Move the boundary to the desired location
Save the edits --> send the data back to the server.

When looking in the PostGIS database for the record and attempting to display the edited polygon it states "unsupported geometries not rendered", previously it would draw the polygon.
Therefore, how do I edit the boundary without causing the geometry error?


texts: Spatial overlay data from many lines into a single line

texts: Hello I have joined lines in my data set into longer lines (based on similar attributes) to reduce the number of rows in my data set. Now I want to take the information from the original smaller lines and update it into the larger lines (two different tables). I would like to concatenate the attributes coming from many lines into the columns for the larger single line. I'm using postGIS for this. thanks


texts: Problem trying to build a stand-alone QGIS application with Unbuntu Trusty system

texts:
I'm currently reading the "Building Mapping Applications with QGIS" Book and I'm stumped when I try to set my paths in Ubuntu.

The first command works fine:
export PYTHONPATH=~/usr/share/qgis/python/

The command to set my QGIS libraries is an issue:
export LD_LIBRARY_PATH=~/qgis path/lib

The issue is that I can't find the "lib" directory.

From what I can see qgis executes from the "/usr/bin/qgis" folder
It's files are located in the "~/usr/share/qgis" folders but there is no hint of a lib folder

I do have the following under the qgis structure:
-doc
-globe
-grass
-i18n
-images
-python
-resources
-svg

Can someone figure out how to set the qgis/lib path under ubuntu trusty?


texts: How to simply report wrong language translation in QGIS?

texts:
I noticed that the "Eliminate sliver polygons" option in QGIS 2.4. has wrong translation in Croatian language (translated as silver), so I decided to report it.

In official page there are instructions how to become translator: http://www.qgis.org/en/site/getinvolved/translate.html#translators-mailing-list
Registration is required for access to mailinglist, so I did it. After registration, there was next message: "Your subscription request has been received, and will soon be acted upon. Depending on the configuration of this mailing list, your subscription request may have to be first confirmed by you via email, or approved by the list moderator. If confirmation is required, you will soon get a confirmation email which contains further instructions."

Few days have passed, and there is nothing in my mailbox.
So, the question is: How to SIMPLY report a wrong language translation in QGIS?
I hope that the Croatian translator will see this.


texts:
I am having the same problem discussed at Where is the attribute table in ArcGIS for Desktop when not visible?.

I don't know what I did but my attribute table is now lurking somewhere but I can't find it. I use only one monitor. The previous question was resolved by right clicking on the taskbar and selecting "show windows stacked." However, I have tried right-clicking everywhere and I haven't seen an option for "show windows stacked."

I've tried quitting ArcMap but nothing changed. I do not want to have to reinstall ArcGIS Desktop nor do I want to have to change the registry.


texts: How to dissolve polygons with common boundary (adjacent) in an Oracle layer through Oracle SQL select command?

texts:
I already have an Oracle layer with polygons which I want to dissolve it through a SELECT SQL command and not use any GIS software, just the Oracle SDO commands.
I can do it normally through GIS tools but I want to do it directly on the SELECT SQL statement.
From this:

To this:

From 3 to 4 polygons, because as you can see I don't want to dissolve polygons that intersects only at a point, they should have a common boundary and not just on point intersected.
I already tried

SDO_AGGR_UNION(SDOAGGRTYPE(geom, 0.005))

but it dissolve all the polygons together at one.


texts: queried layer is not updating after drawing new feature in original shapefile

texts:
I have a project in qgis 1.8 where I added the same shapefile 3 times.  On 2 of them I made a query.  The third one is for edditing purposes.

When I change the value of one of the features (in the third shapefile) it is immediately rendered in the corresponding queried shapefile.  But, when I draw a new feature, nothing changes at all. But when I close and reopen my project, the new feature is drawn in the corresponding queried layer.

Is there a way to update a queried layer with new drawn features, without restarting your project?

thanks


texts: Using MapBox mobile client with GeoServer Protocols?

texts:
MapBox client is beauty and easy to use.

But it seems only support their own protocols (GeoJSON & KML & TileJSON).

GeoServer Only support WCS & WFS & WMS & TMS & WMS-C & WMTS

I am using GeoServer as a map server.

How can I make GeoServer and MapBox works together?


texts: Mapping (Reverb) MODIS NDVI Data in Counties Polygons Shapefiles

texts:
I am quite new to GIS and I have downloaded header files from the reverb website.

My main goal is to get the average NDVI based  on a counties polygon shapefile.


texts: Preserving attributes during kml2shp conversion in ArcGIS for Desktop?

texts:
Looked online and found the two links below on preserving attributes when doing a kml2shp conversion.

Thinking there might be a newer better way of accomplishing this using ArcGIS 10.1 for Desktop?

Similar question asked here 3 years ago: Convert KML to shapefile without losing attribute data

If you convert a complex kml/kmz in ArcMap using ET Geowizard's "Import from Google Earth" tool (or any other kml2shp tool) you will notice many of the informative attributes displayed in Google Earth no longer exist in the resultant Feature Class

Screen shot example

another way to preserve attributes during Convert Kml to Shp


texts: Python package for processing GRIB data in windows

texts:
I would like to know what python packages can be used to process GRIB (Gridded Binary) data on a Windows system. I know there is something called pygrib, but to use that I have to use cygwin windows installation (something I want to avoid).

I am looking for the package to streamline the processing of gridded precipitation data available from NCDC website (for example: http://data.eol.ucar.edu/codiac/dss/id=21.093). NCDC delivers the archived gridded precipitation data as: Unix Compressed (.Z) Gridded Binary (GRIB) Format

Currently I use NOAA Climate and Weather Toolkit (http://www.ncdc.noaa.gov/wct/), but this adds another level of processing (something I want to avoid).

I am seeking python packages, and workflow, perhaps with some examples to automate the seamless processing of GRIB files to derive the precipitation data.