texts: QGIS Upgrading from Madeira (3.4) to Coruna (3.10)

texts:
Quick question!
I am running Madeira, which is my first and only installation of QGIS. Now that
Version Coruna is available.....

Should I bother to upgrade?
If upgrading, should I first Uninstall my
    QGIS Madeira? Or, does the installation routine of the Coruna version take
    care of that (e.g., only over-writing the files that have been
    altered from the one version to another)?
Will the new version of QGIS (Coruna) automatically find my Project files,
and will it carry-over customized projections (coordinate systems)
and so forth?


texts: Extracting coordinates from FeatureCollection to add to Dictionary as values

texts:
I am trying to recreate this ui.Select function with data from an imported FeatureCollection.

In the original example the ui.Select function uses a dictionary that looks like this:

var places = {
  MTV: [-122.0849, 37.3887],
  PEK: [116.4056, 39.9097],
  ZRH: [8.536, 47.376]
};

So I would like to recreate with my own FeatureCollection. My attempt at code for creating the dictionary.

function get_coords(f){return f.centroid()}

var urban_areas = ee.Dictionary.fromLists(boroughs.aggregate_array('lad17nm'),
                                          ee.List(boroughs.map(get_coords)));
print (urban_areas)

My code is here.


texts: Draw the interested point on the map and save the point attribute to OpenGIS database

texts:
I am looking for a map solution can help me pinpoint the wireless route location on my organization's indoor map. I have an indoor map server on which Geoserver is running.

I have found an example page http://openlayers.org/dev/examples/editingtoolbar.html which allows you to draw a point on the map. I also need to edit this point and add attributes for that point such as "Wifi Router Label", "Mac address" etc and save these attributes to PostGIS database so I can have my wireless route geo-location data.

I am lack of javascript knowledge and would like to know which company or person I can contact to get this task done? Or you can post your suggestions here and your help will be appreciated.

Thank you.


texts: What's the purpose of wms_metadata.xml and admin.sld?

texts: I've installed QGIS Server and it works fine. But in /usr/lib/cgi-bin/ folder I've found 2 additional files - wms_metadata.xml and admin.sld. What is the purpose of these file? I've never used them.


texts: Looking for a trajectory similarity measure

texts:
For my final academic work I have developed a compression algorithm for GPS trajectories. I can estimate the quality of the spatio-temporal compression by computing the syncronized euclidean distance (SED) between the compressed and the original trajectory and evaluate my algorithm's performance against a well-known compression algorithm.

A spatio-temporal algorithm, like mine, reduces the trajectory trying to maintain as much temporal information as possible. Spatial algorithms (e.g. Douglas-Peucker algorithm) realize the compression referring only to the spatial characteristics.

What happens now? Considering a spatio-temporal aspect, my algorithm is better than DP. I can assure this by SED measurements. If I plot the three trajectories (original, mine and DP compressed), the trajectory compressed with DP has a better fitting with the original trajectory. The eye-only measurements don't satisfy my need: I need, indeed, an error metric that numerically demonstrates how the DP algorithm is better than mine in a spatial manner.

So I could write:
"Referring to the spatio-temporal factor, my algorithm is better than DP, because it has an SED factor less than DP's SED factor. Alas, simple spatial factor awards DP algorithm because its (name of the new metric) is better than mine".

I have been thinking of perpendicular euclidean distance, but I really don't know if this could be useful. Dynamic Time Warping? What metric could I use for this?


texts: Find csv lat and long points in a shapefile polygon with geopandas spatial index

texts:
I have a csv with two sets of latitude and longitude points, and a shapefile with espg 4326. I want to find which points are inside the shapefile. I've been following the guide here which uses a spatial index to speed up the search however I'm running into an error.

The csv looks like:

e6ed541a288f13745d2755a79372c0c4,fc707fee61baa13b40278a70c679b9a7,3fe94a002317b5f9259f82690aeea4cd,1,2016-11-01T02:03:18.000Z,2,2016-11-01T03:01:13.000Z,2,51.456,0.146,51.506,-0.125,false,false,EE,2,3,2,,,,
36e6a9d201be3957fb78be4dfaacb932,d3cbd057ca3bf69f6805d054b4a228c5,52720e003547c70561bf5e03b95aa99f,1,2016-11-01T16:21:45.467Z,2,2016-11-01T16:43:01.627Z,2,51.499,-0.008,51.506,-0.101,false,false,EE,1,1,1,,,,
f15a4077cc79217e620270a24c127ed4,581c63e4232ae5208f3ae23a09fff3ea,13f9896df61279c928f19721878fac41,1,2016-11-01T20:31:36.000Z,2,2016-11-01T21:21:10.000Z,2,51.529,-0.124,51.512,-0.092,false,false,EI,2,2,1,,,,
bcd71b06e56e2e023fa055a40ee03f20,0ed657e42b4d820eed2affd146067623,3fe94a002317b5f9259f82690aeea4cd,1,2016-11-01T22:58:42.000Z,2,2016-11-01T23:27:24.000Z,2,51.503,-0.113,51.505,-0.086,false,false,EE,2,3,2,,,,

My code so far is:

trip = pd.read_csv(r'test\TripRecordsReporttrips.csv', sep=',',header=None, usecols=[0, 4, 8, 9, 10, 11],names=['TripID', 'Date', 'StartLat', 'StartLon','EndLat','EndLon'])
geometry = [Point(xy) for xy in zip(trip['StartLon'], trip['StartLat'])]
geometry2 = [Point(xy) for xy in zip(trip['EndLon'], trip['EndLat'])]
trip = trip.drop(['StartLon', 'StartLat','EndLon','EndLat'], axis=1)
starts = GeoDataFrame(trip, crs=crs, geometry=geometry) #turn dataframe into geodataframe
ends = GeoDataFrame(trip, crs=crs, geometry=geometry2)
FRC1  = geopandas.GeoDataFrame.from_file('FRC1Shapefile.shp') # import FRC1 polygon

spatial_index = starts.sindex
possible_matches_index = list(spatial_index.intersection(FRC1.bounds))
possible_matches = starts.iloc[possible_matches_index]
precise_matches = possible_matches[possible_matches.intersects(FRC1)]

the error I'm getting is from RTree:

RTreeError: Coordinates must be in the form (minx, miny, maxx, maxy) or (x, y) for 2D indexes

EDIT:
Here's an image of the points sjoin is selecting.


texts: Removing longest segment of linestring and rejoining segments

texts:
I am trying to remove the segment that has the longest length from a linestring and join the two group of segments that remain by the shortest segment between them.
That was easy, using ST_Difference (ST_Split works too) between the linestring and its longest segment, and then I used ST_ShorestLine to join both resulting group of segments.

But I found the biggest problem when the major segment crosses many segments of the linestring. the result of ST_Difference or ST_Split returns more than two groups of segments (because point-type intersections are also taken into account).
This makes the subsequent union of the groups difficult.

What could help me, I think, is that the  ST_Difference or ST_Split (or another function, IDK) only takes into account the intersections of the linestring type but I cannot create the correct query.
Or maybe there is an easier way to remove the longest segment and get the two adjacent linestring to rejoin?
Any ideas?


texts: Using GeoServer WFS to GetFeature by clicking on the map

texts:
Can someone point me to document of crafting a GeoServer WFS GetFeature request with CQL_Filters?

I've already read the GeoServer CQL_Filter Geometric filter section

I have a GeoServer WMS set up to return multiple polygon geometries (from a PostgreSQL database) on the map. I now want to click on the map and in its click event, send a WFS GetFeature request to find the attributes associated with that geometry.

I read the tutorial and it looks like CQL_Filter Contains is what I need but I'm having trouble finding the correct syntax. I'm working with this URL at the moment:

var requestUrl = 'http://localhost:8080/geoserver/wfs?&OUTPUTFORMAT=application/json&REQUEST=GetFeature&SERVICE=WFS&VERSION=1.1.1&srsName=EPSG:4326&TYPENAME=namespace:layer&CQL_FILTER=CONTAINS(geom, Point(' + location.longitude + ',' + location.latitude + '))';

When I post the above URL I get the following error:

<ows:ExceptionReport xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="2.0.0" xsi:schemaLocation="http://www.opengis.net/ows/1.1 http://localhost:8080/geoserver/schemas/ows/1.1.0/owsAll.xsd">
  <ows:Exception exceptionCode="NoApplicableCode">
    <ows:ExceptionText>Could not parse CQL filter list.
Encountered &amp;quot;,&amp;quot; at line 1, column 40.&#13;
Was expecting one of:&#13;
    &amp;quot;-&amp;quot; ...&#13;
    &amp;lt;INTEGER_LITERAL&amp;gt; ...&#13;
    &amp;lt;FLOATING_LITERAL&amp;gt; ...&#13;
     Parsing : CONTAINS(geom, Point(-84.39313136004358,33.757086728151805)).</ows:ExceptionText>
  </ows:Exception>
</ows:ExceptionReport>


texts:
I am trying to show some information on Google Maps which would be in the form of translucent circle overlays. The problem is that some circles overlap, which results in increased opacity and blending of colors if circles are of different color. I searched a lot and came to conclusion that there is no possible way to avoid color blending.

The alternate solution I was recommended was to generete png tiles on my server and display them or decide color of each pixel on server side and use HTML5 canvas and render each pixel. I searched a lot but couldn't find any way to do so. It would be very helpful if someone could point me in the right direction.

I am rephrasing my question below to be more clear and concise

I have to display circles on a map. Radius and color of these circles depict their area of influence and type of influence. These circles are translucent but overlapping should result in only the top one being visible without any blending. How do I generate tiles for this information?

Probably I am not clear so I am adding further details.

Check this jsFiddle  http://jsfiddle.net/Dv4UT/ .There are two overlapping circles green(zIndex : 10) and red(zIndex : 11) with opacity .5 .The overlapping part is a mixture of red and green while I want the overlapping part to be of top color with same opacity ie:red.If this is possible using map overlays then it would be great.Else I think it is possible using HTML5 canvas but that won't work on older version of IE.So finally I would have to fall back on generating tiles on server side and using them as overlays.I don't have any idea of generating tiles and it would be very helpful, if someone could point me in the right direction.I am attaching a screenshot from trulia commute map, check that circles beneath are being overwritten by top ones.I need similar kind of visualization.


texts: Harmonize geometry of different shapefiles (QGIS)

texts:
I have got to harmonize the geometry of a big amount of shapefiles in QGIS

The files are actually very big. So it would be a lot of work to do it manually. Do you have any suggestions for me?


texts: Error while doing Watershed analysis using Python

texts:
I am trying to do watershed analysis using Python in jupyter notebook. Whenever I am passing commands regarding watershed delineation in Python I am getting output as raster matrix and when I am plotting it in graph using imshow() command I am getting squared images in graph and not the desired output.
I am passing commands such as:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import geopandas as gpd
from pysheds.grid import Grid
import mplleaflet
%matplotlib inline
import tifffile as tiff

grid = Grid.from_raster(r'D:\...\ASTGTMV003_N19E073_dem.tif', data_name='dem')
grid.view('dem')

output :
Raster([[202., 204., 204., ..., 517., 524., 532.],
        [204., 202., 202., ..., 514., 519., 523.],
        [202., 201., 200., ..., 519., 516., 519.],
        ...,
        [  0.,   0.,   0., ..., 606., 606., 604.],
        [  0.,   0.,   0., ..., 607., 607., 603.],
        [  0.,   0.,   0., ..., 607., 612., 610.]], dtype=float32)

...
And then while calculating catchment area:
Flow dir:
grid.flowdir(data='inflated_dem', out_name='dir', dirmap=dirmap)
grid.view('dir')

Output:
Raster([[  0.,   0.,   0., ...,   0.,   0.,   0.],
        [  0.,   2.,   2., ...,  -1.,   4.,   0.],
        [  0.,   4.,   1., ...,  64.,  -1.,   0.],
        ...,
        [  0.,   4.,   4., ...,   1., 128.,   0.],
        [  0.,  -1.,  -1., ..., 128.,   1.,   0.],
        [  0.,   0.,   0., ...,   0.,   0.,   0.]], dtype=float32)

But then I plotted it in image using:
plotFigure(grid.dir,'Flow Direction','viridis')

And I got image but while calculating catchment area
By specifying x and y:
x, y = 73.6, 19.6
grid.catchment(data='dir', x=x, y=y, dirmap=dirmap, out_name='catch',
               recursionlimit=15000, xytype='label')
grid.clip_to('catch')
catch = grid.view('catch')
plt.imshow(catch)

I am getting output as

Why is this occurring?


texts: Avoiding arcpy.mapping.ExportToPDF locking mxd with data driven pages?

texts:
I've written a script that makes use of a data driven pages mxd, but I'm having issues with locks. The script iterates through a series of feature classes and uses them as the base layer for the DDP. I've found that once an mxd object is created, saved, and the code arcpy.mapping.ExportToPDF is used, the mxd can't be saved again. The code mxd.save () fails on its second iteration. I've run the script commenting out arcpy.mapping.ExportToPDF and the script ran fine. Here's my code:

oldMapDoc = mapDoc

#iterate feature classes
for i, fc in enumerate (sorted (fcs)):
    print fc
    mxd = arcpy.mapping.MapDocument (mapDoc)
    df = arcpy.mapping.ListDataFrames(mxd)[0]
    ddpLyr = [l for l in arcpy.mapping.ListLayers (mxd) if l.name == "Symbology"] [0]
    ddpLyr.replaceDataSource (platsGdb, "FILEGDB_WORKSPACE", fc)
    mapFolder = os.path.join (distFolder, fc.title ().replace ("_", " "))
    mxd.save ()
    pageCount =  mxd.dataDrivenPages.pageCount
    print pageCount
    for pageNum in range(1, pageCount + 1):
        mxd.dataDrivenPages.currentPageID = pageNum
        index = mxd.dataDrivenPages.pageNameField.name
        plat = mxd.dataDrivenPages.pageRow.getValue (index).replace ("\\", "_")
        outPdfName = fc + "_" + plat + ".pdf"
        pageId = mxd.dataDrivenPages.currentPageID
        print "Exporting page {} of {}: {}".format(pageId,
                                                   pageCount,
                                                   outPdfName)
        outPdf = os.path.join (mapFolder, outPdfName)
        if arcpy.Exists (outPdf):
            print "already made"
            continue
        arcpy.mapping.ExportToPDF (mxd, outPdf)
        #break
    del mxd, df, ddpLyr, index, plat, outPdfName, pageId, mapFolder, pageCount
    mapDoc = oldMapDoc.replace (".mxd", "") + "_" + str (i) + ".mxd"
    print "copying"
    copyfile (oldMapDoc, mapDoc)

I've tried the variation of creating the mxd object outside of the loop. This produces a different error - the script runs but the DDP layer never actually updates.

As you can see, I've created a workaround where I copy the mxd file with each iteration. It works, but it's pretty wonky. Has anyone dealt with this issue before and found a solution besides something hokey like I have?


texts: openlayers3: layerswitcher + getfeatureinfo

texts:
I'm working on a webmapping project which use geoserver and OpenLayers 3 for the vizualisation. I'm stuck when I try to use "WMS GetFeatureInfo" and "LayerSwitcher". Both functions work individually but when I try to combine both it doesn't work anymore. I'm sure I'm on the right way but missing something...

For the GetFeatureInfo I use code like in this example: http://openlayers.org/en/v3.0.0/examples/getfeatureinfo-tile.html

and for the LayerSwitcher I have implemented this: https://github.com/walkermatt/ol3-layerswitcher

Do I may have to change something in the ol3-layerswitcher.js ?

    <link rel="stylesheet" href="web/css/ol.css" type="text/css">

    <link rel="stylesheet" href="web/css/ol3-layerswitcher.css" />

      <title>test</title>

&nbsp;

At the end I like to have a Layerswitcher which can group my layers and it should be possible to use "GetFeatureInfo" for every Layer.


texts: QGIS 3D - no terrain altitude binding for polygons / buildings

texts:
My first post to this topic:
QGIS 3D - no terrain altitude binding
There I have had problems with 3D objects from a point layer. I solved this.

Now I have a problem with 3D objects from a polygon layer.
I use polygons / buildings from osm. All buildings with the same high. Within qgis2threejs is all ok. But with QGIS 3D not: No altitude binding.
But I need for my ideas 3D of QGIS, I can not use qgis2threejs.
Any ideas to solve the problem? I need the same result as for the buildings within qgis2threejs.

I get only a good result if I use flat terrain.

Here a screenshots of the problem and of my settings:

above: The buildings with the black lines are examples, not the only wrong building-heigh-positions.
The buildings left are OK, on the right wrong.

My settings:

Without terrain: no problems

Windows 10, QGIS 3.10.6

edit: try on QGIS 3.12.3: only crashes....oh....


texts: Avoiding Error 1325 EsriPyFldr not valid short file name when uninstalling ArcGIS for Desktop?

texts:
I am trying to uninstall ArcGIS 10.3.1 for Desktop from Windows 7 SP1 using Uninstall a program and I am receiving an error message (pictured below):

  Error 1325. EsriPyFldr is not a valid short file name.

The error message dialog appears just after:

  Computing space requirements

seems to complete.

If anyone else has encountered this, then have you been able to find a way to get past it and uninstall ArcGIS 10.3.1 for Desktop?

My machine has what I would call a vanilla install of ArcGIS 10.3.1 for Desktop and ArcGIS Pro 1.1.1 so there are Python versions 2.7.8 and 3.4.1 as a result. ArcGIS Earth Beta is also on the machine.

I have made a posting on the same issue to https://geonet.esri.com/message/575655#575655

After running this utility (https://support.microsoft.com/en-us/mats/program_install_and_uninstall) unsuccessfully (it hung for two hours before I killed it) the error message when trying to uninstall shifted to:

  This action is only valid for products that are currently installed

ArcMap 10.3.1 continues to start fine and doing
wmic product where "Vendor like '%Environmental Systems Research Institute%'" get Name, Version, InstallDate, InstallLocation in a DOS window (at the suggestion of @MattWilkie) reports only the presence of Python 3.4.1 for ArcGIS Pro, ArcGIS Pro and ArcGIS Earth on my machine.


texts: How to Avoid Object ID conflicts in Versioned geodatabase

texts:
We have just implemented a versioned Geodatabase and our schema has few feature classes and relationship classes.

We do not Reconcile and Post to the parent ( Default) version, as we need some intermediate monthly job extraction and do not want to mix up jobs of past month.

Our concern is , since we do not Post to default on regular base, is there a chance that Object Ids be duplicated in the monthly versions as the same as Default ? In that case wont we lose features once we post to the default at the end of the month ?

Or will the Oracle/ArcSDe take care of the Object ids for all versions , so that the basetable( Default) and the delta tables( Versions) does not have features with same Object ID in a feature class ?


texts: Processing array of values at cell location across raster maps using ArcGIS Desktop

texts:
I have several (say n=10) raster maps with ratings values. At a particular cell location, I would like to create a sum composed of the maximum rating, half the next highest rating, then a quarter of the next highest, etc.

I understand how to find the MAX, MIN, AVERAGE ... across maps, but not how to order the values at a cell location to process them.

I am using ArcGIS Desktop.

Suggestions?


texts: Locate LAS Points By Proximity gives ERROR 999999?

texts:
I am trying to locate LAS points which is within a specific distance (25cm or 50cm) from the 3D points in a feature class. This is exactly what the Locate LAS Points By Proximity tool should do. My problem is just that I get a 999999-error, and I have no clue why.

My input files are a Las Dataset and a 3D point feature class and I have tried a lot of combinations of settings and outputs, but with the same error message every time.

I can only drag-and-drop my point feature class into the tool - if I browse to find it in the windows folder, then it is not visible. Maybe that can be a hint for solving the issue.


texts:
I have a raster file, used "clip raster by mask layer", then I used rastercalculator to set the values of the new rasterlayer to zero and now I would like to merge them together again.

However, the merge tool doesn't work and it doesn't give any error messages except telling me that it's not working.

The raster files have the same projection, both are GEOTIFF, they have the same resolution and both are Float 32... perhaps there is something I'm not seeing?

Here I include images of the used settings and the output dialog:

The red error message just says "the following layer has not been created. More information on how this algorithm works can be found in the QGIS main Window".

I'm using the standalone LTR Version QGIS 3.4.14

What could be a reason why it's not working?


texts: Finding model created using ArcGIS ModelBuilder?

texts:
I made a model in the ModelBuilder (arcgis 10.2) and saved it. Unfortunately, not some place where I can find it again.

However, if I want to save a new model (save as..),the old model is there in toolbox.tbx.

Does anyone know how I can re-open my model again? ArcCatalog is not helping me so far.


texts: Understanding Merge Error 001156 from ArcGIS Pro?

texts: I have been having trouble merging point features together in ArcGIS Pro due to the error "001156." I've been looking up solutions but they don't necessarily help me out! I don't know which output value to use either...ArcGIS Pro is fairly new to me!


texts: ArcMap data links redirecting to location of exported map

texts:
This issue had only started in recent days. Just wondering if anyone else has come across this before and might know of a solution?

I have multiple mxds that contain layers that normally live in a folder called "Data" within a project folder.

So the source for these shapefiles would look something like:
D:\Project\Data

The mxd would live here:
D:\Project

Then there would be a separate folder for where the maps (JPEGs/PDFs etc.) are exported out to:
D:\Project\Maps

What seems to be happening with some of the mxds, is that the links for the data are redirecting to the maps folder after you've closed the mxd and reopened it the next day.
So instead of going to D:\Project\Data, when you reopen the mxd, the data is linking to D:\Project\Maps\Data.

We are quite certain that this is only occurring when we've just exported a map and saved the mxd (whether the order is save then export, or export then save, we're not 100% on).

This has been happening at multiple workstations with project folders and data located in different locations.
All mxds have "Store relative pathnames to data sources" checked as well.


texts: I am importing points containing elevation from a survey controller into a GIS. The survey template on the controller contains a projection using a HT2 Geoid model, in NAD 83 Canada UTM. I can match the EPSG to 26914 in a GIS, which has the same spheroid parameters but recently I have noticed that the flattening is different between the template in the controller and EPSG:26914 in GIS. Can anyone tell me what impact the difference in flattening will have? My guess is, because the elevation is recorded on the controller it won't change regardless of what geoid model parameters are in place in the GIS.


texts: Overwrite previously exported Excel file

texts:
I am exporting a table to Excel and I would like the overwrite access to be on. Below is the code line. Please let me know how I can give overwrite access through the code if it is possible.

#Export to Excel
arcpy.TableToExcel_conversion("UniqueTandR",toolpath+"/FinalResult.xls","ALIAS","DESCRIPTION")


texts: Importing GeoPackage via GeoPackage vs. Import Vector in QGIS

texts:
Is there any difference when importing a GeoPackage into QGIS using the data source manager between using the section "Vector" and the section "GeoPackage"?
I did not see any difference so far, but this may be because my geopackages are rather simple.


texts: Table to Table works differently when implemented in Jupiter Notebook and ArcMap

texts:
I tried to import a lot of tables from geodatabase to .csv to later work with them. I encountered strange values and decided to do it manually for one table through the ArcMap interface. Everything went alright there. I then copied the tool as a Python snippet, ran it in Jupiter and got the same (wrong) values as the first time.

To clarify: I run the same code in ArcMap and Jupiter and get different results

arcpy.TableToTable_conversion(in_rows="mean_1312", out_path="C:/.../Pre And Tmp", out_name="tmp_mean_1312.csv")

Strangely, not all columns are different between these two versions.

Result from Arcmap:

Result from Jupiter:


texts:
I have vector created by conversion of a raster to vector operation. I want to clip this vector file.
While I am trying to do that, I am getting the below mentioned error:
Feature (355) has invalid geometry. Please fix the geometry or change the Processing setting to the "Ignore invalid input features" option.
Execution failed after 0.05 seconds
What am I doing wrong?


texts: Implementing IEditSketch interface for custom CAD tools?

texts:
I want to build cad tools which will be used in municipality applications on esri.

The user wants to edit on the map and does it by using IEditSketch.But i want to write a custom IEditSketch.So i want to implement that interface.Can u share your experiments about it.


texts: SpatialPolygonsDataFrame chloroplethr not writing to KML for GoogleMaps using R

texts:
I'm trying to do a simple chloropleth map for postal codes. I have the orginal kml exported here.

When importing into R, i'm doing:

library(rgdal)
library(rworldmap)
library(plotKML)
library(maptools)

foo <- readOGR("CP_04Camp_v2.kml")

This yields a SpatialPolygonsDataFrame :

class       : SpatialPolygonsDataFrame
features    : 176
extent      : -92.46879, -89.38602, 17.81287, 20.455  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
variables   : 2
names       : Name, Description
min values  :     ,
max values  :     ,

First, I tried to add a column in the foo@data with the color attribute, but I keep having a hard time exporting the kml to use in Google Maps.

I have been trying:

> plotKML::kml(obj = foo, "test.kml")
KML file opened for writing...
Error in FUN(X[[i]], ...) :
  Proj4string does not contain \+datum= parameter.
 Consider converting to the referent CRS +proj=longlat +datum=WGS84 manually.

Or using maptools:

> kmlPolygons(obj = foo, kmlfile = "test.kml")
Warning message:
In rep(col, length.out = length(obj)) :
  'x' is NULL so the result will be NULL

I get the same error messages even if I don't even "touch" the imported object. What am I missing?

session info:

> sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] maptools_0.8-41 plotKML_0.5-6   rworldmap_1.3-6 needs_0.0.3     rgdal_1.2-5     sp_1.2-4        xlsx_0.5.7
 [8] xlsxjars_0.6.1  rJava_0.9-8     readxl_0.1.1    ggmap_2.6.1     RODBC_1.3-14    memisc_0.99.7-1 MASS_7.3-45
[15] lattice_0.20-34 magrittr_1.5    ggplot2_2.2.1   inegiR_1.2.0    tse.utils_0.1.0 tse.sql_0.1.0   stringr_1.1.0
[22] dplyr_0.5.0     foreign_0.8-67

loaded via a namespace (and not attached):
 [1] nlme_3.1-128        pixmap_0.4-11       pbkrtest_0.4-6      xts_0.9-7           devtools_1.12.0
 [6] RColorBrewer_1.1-2  backports_1.0.4     tools_3.3.2         R6_2.2.0            rpart_4.1-10
[11] Hmisc_4.0-2         DBI_0.5-1           lazyeval_0.2.0      mgcv_1.8-15         colorspace_1.2-7
[16] nnet_7.3-12         raster_2.5-8        withr_1.0.2         gridExtra_2.2.1     quantreg_5.29
[21] htmlTable_1.8       gstat_1.1-4         SparseM_1.74        checkmate_1.8.2     scales_0.4.1
[26] classInt_0.1-23     digest_0.6.10       minqa_1.2.4         colorRamps_2.3      RSAGA_0.94-5
[31] dismo_1.1-4         htmltools_0.3.5     base64enc_0.1-3     jpeg_0.1-8          lme4_1.1-12
[36] plotrix_3.6-4       maps_3.1.1          FNN_1.1             zoo_1.7-13          jsonlite_1.2
[41] acepack_1.4.1       car_2.1-4           Formula_1.2-1       shapefiles_0.7      geosphere_1.5-5
[46] Matrix_1.2-7.1      Rcpp_0.12.7         munsell_0.4.3       proto_1.0.0         stringi_1.1.2
[51] plyr_1.8.4          grid_3.3.2          parallel_3.3.2      splines_3.3.2       mapproj_1.2-4
[56] knitr_1.15.1        rjson_0.2.15        spacetime_1.2-0     reshape2_1.4.2      XML_3.98-1.5
[61] latticeExtra_0.6-28 data.table_1.10.0   png_0.1-7           spam_1.4-0          nloptr_1.0.4
[66] aqp_1.10            RgoogleMaps_1.4.1   MatrixModels_0.4-1  gtable_0.2.0        reshape_0.8.6
[71] assertthat_0.1      e1071_1.6-7         class_7.3-14        survival_2.40-1     tibble_1.2
[76] intervals_0.15.1    memoise_1.0.0       fields_8.10         cluster_2.0.5


texts: What signal does x button in messageBar emmit or how to get rid of it using PyQGIS?

texts:
I'm creating QGIS plugin and I'm informing user about the progress of computation via QProgressBar, that is shown in messageBar. I have added Cancel button, that gives user a possibility to end the computation (it shows QMessageBox with question whether the user wants to end computation or not). I'd like to implement the same functionality as Cancel button to x button, that is created automatically as part of messageBar or completely erase it, so there will be only Cancel button.

Does anybody have an idea how to get rid of x button or what signal is emitted, when it is clicked? I found one signal in QObject class - destroyed, but if I understand it correctly from the documentation, the process of destroying QObject cannot be stopped this way.

For illustration, my code looks like this:

def progressBar(self):
    """Initializing progress bar.
        """
    self.progressMessageBar = iface.messageBar().createMessage(u"Ground Radiation Monitoring:",u" Computing...")
    self.progress = QProgressBar()
    self.progress.setMaximum(100)
    self.progress.setAlignment(Qt.AlignLeft|Qt.AlignVCenter)

    self.cancelButton = QtGui.QPushButton()
    self.cancelButton.setText('Cancel')
    self.progressMessageBar.layout().addWidget(self.cancelButton)
    self.progressMessageBar.layout().addWidget(self.progress)
    self.iface.messageBar().pushWidget(self.progressMessageBar, iface.messageBar().INFO)

    self.cancelButton.clicked.connect(self.onCancelButton)

And messageBar:


texts: Error 999999 in Zonal Statistics estimation using ArcPy?

texts:
I have multiple rasters saved in "Zonal.gdb". All the rasters in "Zonal.gdb" and "TownWilson_PrjUTM" have projected coordinates "WGS_1984_UTM_Zone_15N". I want to clip each raster to the extent in the shapefile (TownWilson_PrjUTM) and resample to 10 meter resolution and estimate zonal statistics as shown in the code below.

I am getting an error as shown below. Cannot understand why I am getting the error even though I have projected rasters and shapefile to the same projected coordinates.

Can somebody help me to figure this out?

import arcpy, os, arcinfo
from arcpy import env
from arcpy.sa import *

arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
env.workspace = "H:/Arkansas_Flood/Zonal.gdb"
watershedFeat = "H:/Arkansas_Flood/Zonal.gdb/TownWilson_PrjUTM"
outDir = "H:/Table/"

for raster in arcpy.ListRasters():
    raster_name = os.path.basename(raster).rstrip(os.path.splitext(raster)[1])
    ## Clip Raster Dataset by known extent - Left Bottom Right Top
    clip_raster = arcpy.Clip_management(raster, "231366.658400 4140371.653300 280675.622100 4180783.072200", "#", "#", "#", "NONE") # with TownWilson_Prj
    # Process: Project Raster
    finalras = arcpy.Resample_management(clip_raster, "#", "10", "BILINEAR")
    outTable = outDir + raster_name + "_TBL.dbf"
    arcpy.gp.ZonalStatisticsAsTable(watershedFeat,"T_R", finalras, outTable,"NODATA","MEAN") # DATA for precipitation, NODATA for RZSM

arcpy.CheckInExtension("Spatial")

> Error:
>
> Executing: zonalMeanRev Start Time: Thu Jan 04 23:50:02 2018 Running
> script zonalMeanRev... Failed script zonalMeanRev...
>
>     Traceback (most recent call last):   File "C:\Users\jayaskeradl\Desktop\Arkansas_
> Flood\My_Phyton_Scripts\zonal_mean_rev.py", line 46, in <module>
>     clip_raster = arcpy.Clip_management(raster, "231366.658400 4140371.653300 280675.622100 4180783.072200", "#", "#", "#", "NONE") # with TownWilson_Prj   File "c:\program files
> (x86)\arcgis\desktop10.4\arcpy\arcpy\management.py", line 13771, in
> Clip
>     raise e ExecuteError: ERROR 000622: Failed to execute (Clip). Parameters are not valid. ERROR 000628: Cannot set input into
> parameter in_raster.
>
>
> Failed to execute (zonalMeanRev).
>
>
> Traceback (most recent call last):   File
> "C:\Users\jayaskeradl\Desktop\Arkansas_
> Flood\My_Phyton_Scripts\zonal_mean_rev.py", line 46, in <module>
>     clip_raster = arcpy.Clip_management(raster, "231366.658400 4140371.653300 280675.622100 4180783.072200", "#", "#", "#", "NONE") # with TownWilson_Prj   File "c:\program files
> (x86)\arcgis\desktop10.4\arcpy\arcpy\management.py", line 13771, in
> Clip
>     raise e ExecuteError: ERROR 999999: Error executing function. The table already exists. No spatial reference exists. Failed to execute
> (Clip).
>
>
> Failed to execute (zonalMeanRev).

Traceback (most recent call last):
  File "C:\Users\jayaskeradl\Desktop\Arkansas_ Flood\My_Phyton_Scripts\zonal_mean_rev.py", line 53, in <module>
    finalras = arcpy.Resample_management("H:/Arkansas_Flood/Zonal.gdb/newRaster", "#", "10", "BILINEAR")
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\management.py", line 14243, in Resample
    raise e
ExecuteError: ERROR 999999: Error executing function.
The table already exists.
No spatial reference exists.
The table was not found. [VAT_newRaster]
Failed to execute (Resample).

Failed to execute (zonalMeanRev).


texts: Separating north, south, east and west of each parcel according to their neighborhood using ArcGIS Desktop?

texts: I am searching an easy way, like code, I have about 70000 parcels from the city. How can I separate north, south, east and west of each parcel according to their neighborhood?


texts: Closest point on a line segment from point outside the line

texts:
From the above pic I found the value x from line (p1,p2) and point a.

Using y=mx+b and an imaginary red line which is perpendicular to the black line having slope -1/m

and the intersecting point x, the problem is that x is not actually intersecting.

See pic below


texts: Seeking formula for Distance Calculation on Mars

texts:
I have a database of craters on Mars in decimal latitude and longitude. The latitude ranges from about -90 to 90 and the longitude from -180 to 180 which feels right to me.
Is there a formula to calculate the distance between two points using this system on Mars?
I am using data from: Robbins Crater Database Explanation of the fields can be found at: A new global database of Mars impact craters ≥1 km


texts:
I have a table with three attributes [field1, field2, field3].

I'm trying to insert a new feature to this table by defining field2 and field3 only.

lyr = QgsMapLayerRegistry.instance().mapLayer('layer')

new_feat = QgsFeature()

new_feat.setAttribute(1,"text1")
new_feat.setAttribute(2,"text2")

lyr.dataProvider().addFeatures([new_feat])
lyr.commitChanges()

It seems that I have a problem with the index in the setAttribute.
What am I missing?


texts: Using the pixel quality layers of the Landsat SR products to create cloud and water masks in R

texts:
I tried to use the RStoolbox package function classifyQA to generate information on quality of the layers of interest, namely - cloud, cloud shadows, and water. I am unable to understand how to mask clouds and water using this function.
#One of the attempts
cloudMask(img,crop_pixelqaband_1992) # trying to use the pixel QA band provided by SR products to create a mask. I understand right now the code is not specific in pulling out only the cloud and cloud shadow data from the QA band layer.

Return error
Error in .local(x, ...) : invalid layer names # the function only seems to take blue and thermal bands

May be I am doing this wrong? Is there another way to do cloud and water masking using Landsat surface reflectance products that do not provide thermal bands in R? Can someone please suggest how one can go about carrying out masking for clouds and water in R using the pixel quality band of SR products?


texts:
Please recommend pgrouting tutorials or workshops that are most current.

I ask because I just read a post that claims a tutorial is deprecated. (see osm2po vs osm2pgrouting: workshop post)

I am working through the pgrouting-workshop (http://workshop.pgrouting.org/) and find it mostly helpful.

The OS and environment here:  ubuntu 12.04, postgres 9.1, postgis 1.5.

While it may be unreasonable to ask that all deprecated tutorials be so marked, it would help if the date of authorship is included prominently in the header of all training materials, then newbies would know to expect some staleness.


texts: Splitting overlap between polygons and assign to nearest polygon using R?

texts:
My question is how to perform the process described in the answer to this similar question in R, ideally using the sf package. In the linked case this is done using using PostGIS with some sf commands, but I am looking for a native R solution.

Basically: given a set of overlapping (multi)polygons, how do I split the overlapping polygons according to the non-overlapping sections which are closest.

Here is a reproducible example along with my progress so far (with help from the r-spatial vignettes and github).

library(sf)
library(tidyverse)

pol = st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))))
b = st_sfc(pol, pol + c(.8, .2), pol + c(.2, .8))
par(mar = rep(0, 4))
plot(b, col = NA)

i = st_intersection(st_sf(b))
par(mar = rep(0, 4))
cl = sf.colors(3, categorical = TRUE)
plot(st_geometry(b))
plot(st_geometry(i[i$n.overlaps == 3,2]), col = cl[1], add = TRUE)
plot(st_geometry(i[i$n.overlaps == 2,2]), col = cl[2], add = TRUE)

d = st_difference(b)
plot(d, col = cl)

So far I have been able to create the divisions I want at the voronoi partitions among the centroids of each polygon:

independent <- b %>% st_sf %>% st_intersection %>% subset(n.overlaps<=1)
overlap <- b %>% st_sf %>% st_intersection %>% subset(n.overlaps>1) %>% st_union()
partition <- b %>% st_centroid %>% st_union %>% st_voronoi %>% st_cast %>% st_intersection(overlap)

plot(st_geometry(independent), col=cl)
plot(st_geometry(partition), col=cl, add=TRUE)

How do I assign each partition of the overlap to the appropriate polygon from the non-overlapping sections? Is this the best method?


texts: How does the QGIS construct class intervals?

texts:
I'm using graduated style of a point layer to assign a colour ramp to symbols (squares). I assigned manual values to class intervals specifying classes to be:

0.60 - 0.80
0.80 - 1.00

That is, in entering class bounds for a lower class I used 0.600000 as lower value and 0.800000 as upper. And then 0.800000 as lower for next class.

How are then the class intervals constructed in such case? Are they:

(0.60 - 0.80)
[0.80 - 1.00)

Or:

(0.60 - 0.80]
(0.80 - 1.00]
...


texts:
is it possible to enhance the shadow effect to choose the value from the table instead of mm or pixel? maybe in the next version of QGIS.
Or is there any other workaround


texts: How to prevent popup when there is no feature

texts:
I use OL3 together with geoserver to display one WMS layer. Then I use the GetFeatureInfo to display some of the data on button click. I want to prevent to show the popup when I click outside the layer's features (as shown in the attachment)

How can I do this in OL3?
This is the part of the code, which takes care the getFeatureInfo:

map.on('singleclick', function (evt) {

  var viewResolution = map.getView().getResolution();
  var coord = evt.coordinate;
      if (grid_layer.getVisible() && grid_layer.get('name') != 'Basemap') {

          var url = grid_layer.getSource().getGetFeatureInfoUrl(coord, viewResolution, 'EPSG:3857', {
              'INFO_FORMAT': 'text/html',
              'FEATURE_COUNT': '1',
              'propertyName': 'url'
          });

          if (url) {
              overlay.setPosition(coord);
              content.innerHTML = '<iframe src="' + url + '" style="border:0px #BD8D46 dotted;" scrolling="auto" frameborder="no" align="center"></iframe>';

          }
      }
  });


texts: Unify the color of graduated classes across several maps

texts:
I'm using ArcGIS 10.4 (advanced license).

The two maps below don't have the same minimum and maximum values.

I have tried using Save As layer file > Import symbology but it didn't work?

How is it possible to unify the graduated colors across the two layers?


texts: OpenLayers 'zoom to extent' issue on swipe tool

texts:
I am working with OpenLayers version 3.16 in which I have added few layers and comparing those using swipe tool OL3-ext. Layers are coming from GeoServer.

While zooming to extent on selected layers my OpenLayers 3 is crashing. Please check attached screenshots as following.

Tried flowing to resolve error:-

1.In single map (on one layer) zoom to extent working fine.

2.Updated my layers order.

3.Checking related questions in GIS.SE

Screenshots:-

1.Before comparing two different layers using swipe tool selected layer highlighted in red color.

2.After zoom to selection OpenLayers crashes (highlighted in red color)


texts: Where does Geo-Platform fit in a GeoStack?

texts:
I am totally new to GIS development and I am trying to grasp the concept.

From my understanding the stack looks as follows (Please correct me):

Back End GIS Server: GeoServer/GeoNode/ArcGIS Server (Some of the back end servers do have their own spatial database - however, PostGIS can still be used in conjunction with them)

Back End Spatial Database (optional): PostgreSQL/PostGIS database can be used in conjunction with back end GIS server (ArcGIS for example)

Front End Library (to display the shapes/points on top of a map(Google Maps/Open Street/etc.)): Openlayers, lefletjs and GeoExt

Now I am reading about Geo-Platform and I am not sure where exactly does it fit in the above?


texts:
I am working with a GIS problem using a single input of a polygon shapefile.

Consider an irregular polygon. I want to draw vertical lines across the extent of the polygon at equal spacing.

How I intend to proceed is:

Identify the bounding box (done using PyShp)
Draw vertical Lines parallel to the left edge of the bounding box at equal spacing (How?)
Clip the lines to the extent of the polygon (How, without using ArcPy?)

Note: They are required to be only vertical, and not a graticule. Also, I do not intend to use ArcPy, and intend to perform the coding completing in Python (2.7) as this segment of code needs to go into a tool generated from PyQt.


texts:
I have point data in GDA2020.  shapefile *.prj file says units are meters.  when I load the data in QGIS the layer units are degrees nd the buffer tool is always degrees.  I have searched for every setting I can find yet the layer and buffer tool are always degrees.  I have exported and re-imported, restarted QGIS.  regardless of the CRS settings the map units are always degrees.  What am i missing ?

Source data came from a table containing lat/long it looks fine in tableau or even in QGIS.  I just can't buffer it because it's not meters.  Need to do a spatial select on another point feature set.  Also, data set is Australia wide


texts: I am trying to find a way to locate any unsnapped pipes/manholes. I noticed that my snapping environment was not set up in ArcMap the first few times I made edits and I want to go back and reverse those mistakes. There are thousands of pipes and manholes so I am sure I am not the only one to have this happen.


texts:
I would like to use the Positron base map provided by CartoDB with CesiumJS.
According to this blog post this should be possible using this cesium-cartodb plugin. I'm using Cesium version 1.6.

Though the example given on the plugin github README page appears very straight forward, it is not clear what URL to use. I can't find what the values {s}, {z}, {x} and {y} relate to. I would also like to know if this image provider will automatically load in appropriate zoom level tiles like the default CesiumJS image providers.

var basemapProvider = new Cesium.CartoDBImageryProvider({
    url: '<MAP_TEMPLATE_URL>', // e.g.:     http://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png,
    credit: 'Basemap courtesy of CartoDB'
});

viewer = new Cesium.Viewer('cesiumContainer',{
    imageryProvider: basemapProvider,
    baseLayerPicker: false
});

None of the examples ive seem work for me, they make calls to the following urls and fail:

http://b.basemaps.cartocdn.com/light_all/1/1/0.png/
http://b.basemaps.cartocdn.com/light_all/1/0/1.png/
http://c.basemaps.cartocdn.com/light_all/1/1/1.png/
http://a.basemaps.cartocdn.com/light_all/1/0/0.png/


texts: Display millions of locations very, very fast

texts:
We're building a mapping application for voter data using PostGIS and Leaflet. It has a web interface. Our voter file will have several million entries. In our app you can select a subset of the voter file and display it.

It works fine when you're zoomed in because the result set is small. The challenge is to figure out what to do when the result set is large. As I see it, we have multiple options:

Just bail out when result sets exceed X records. The user would receive an error message: "Too many records to display. Zoom in or narrow your search criteria." This isn't ideal.
Try to do client-side clustering and just display the clusters. This won't work because it's not possible to transfer very large result sets to the browser quickly. No amount of compression will get a million records down to an acceptable size.
Try to do server-side clustering. We've done that using the techniques found here: Spatial clustering with PostGIS? The difficulty is, again, speed; you can't aggregate millions of locations and calculate cluster centroids and counts quickly. Maybe we're doing it wrong.
Do some kind of aggregation using zip codes. This is faster; essentially you do a select zipcode, count(*) from voters where (some predicate) group by zipcode. You then look up the zipcode centroid separately. It's faster because it's a covered query which can usually be answered by hitting indexes alone. But zipcodes don't display well because they're dense in a city and not so dense elsewhere.
Aggregate based on some other field, maybe something related to tiles. I haven't fully thought this one through.
Exploit some property of the GiST index on the location. This is theoretically possible; a GiST index already "clusters" points internally in the R Tree. I don't know where to start with this idea.

How do other people solve the problem of summarizing millions of points, dynamically, quickly, for presentation to an end user?


texts:
I know similar questions have been asked several times, but none of the solutions seems suitable for my problem or don't work..

I computed flowlines using the r.flow tool from GRASS GIS, then I export these using v.out.ogr which gives me the files flowlines.dbf, flowlines.prj, flowlines.shp and flowlines.shx. The image shows the figure in GRASS of my flowlines over a DEM.

Now I need to read this shapefile into python and reproduce that same image in python. To read the file I downloaded the shapefile module, which seems to work fine, but after that I'm lost. I was told to use fiona module by a friend but it doesnt want to install on my ubuntu machine, not sure what's the problem, so a solution without fiona would be great..

Does anyone have good advice?


texts: How to convert a .csv file to shapefile with the attribute table geometry using pyqgis

texts:
I am trying to convert a .csv file to .shp using pyqgis. With this publication ETA managed to generate the file with all its extensions and the table of attributes but I do not see the entities (type point). column 1 and 2 have the information of Lat and Long the code I have is this:

ruta = "C:/Users/Pc/Desktop/MANUELITA/productividad/CSV/35300030000_C58_b.csv" # set the filepath for the input CSV
inp_tab = QgsVectorLayer(ruta, 'Input_Table', 'ogr')
salida = 'C:/Users/Pc/Desktop/MANUELITA/productividad/CSV/prueba.shp'
fields = inp_tab.fields()
crs = 4326
outLayer=QgsVectorFileWriter.writeAsVectorFormat(inp_tab, salida,"fields" ,inp_tab.crs() ,"ESRI Shapefile", layerOptions=['SHPT=POINT'])

the .csv fields looks like this:


texts: Getting historical imagery of Google Maps?

texts: I want to show historical images of google map of particular location and images will change when we click previous next button according to its published date ?


texts: Creating grid in QGIS representing trial plots

texts:
I would like to create a grid in QGIS. This grid shall represent agricultural trial plots (each plot about 4,5 x 10 m). Digitizing each polygon is time-consuming and not very exact.
Is there a simple way to do that?
I tried "create grid" but that always aligns the grid exactly north/south and that is not what I need. I would like to place the grid off that direction.
There was a plugin "breeder map" but it does not seem to be available (QGIS 3.10). Any ideas?


texts:
I am quite new to GIS and qGIS especially.

I would like to install cfarmers cartogram plugin. From what I understand, this has to be done manually - so I 've followed the instructions from RobinLoveLace (see link below).

Where has the cartogram plugin for QGIS gone?

When I try to click the icon for the Cartogram Creator, I get an error message (see below).

Any tips on how to get it to work?

  An error has occured while executing Python code:

  Traceback (most recent call last):
    File "C:/PROGRA~2/QGISCH~1/apps/qgis/./python/plugins\cartogram\plugin.py", line 136, in run
      d = doCartogram.Dialog(self.iface)
    File "C:/PROGRA~2/QGISCH~1/apps/qgis/./python/plugins\cartogram\doCartogram.py", line 61, in init
      for name, layer in layermap.iteritems():
    File "C:/PROGRA~2/QGISCH~1/apps/qgis/./python/plugins\cartogram\doCartogram.py", line 61, in init
      for name, layer in layermap.iteritems():
    File "C:\PROGRA~2\QGISCH~1\apps\Python27\lib\bdb.py", line 49, in trace_dispatch
      return self.dispatch_line(frame)
    File "C:\PROGRA~2\QGISCH~1\apps\Python27\lib\bdb.py", line 68, in dispatch_line
      if self.quitting: raise BdbQuit
  BdbQuit

  Python version:
  2.7.4 (default, Apr  6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)]

  QGIS version:
  2.4.0-Chugiak Chugiak, 8fdd08a

  Python path: ['C:/PROGRA~2/QGISCH~1/apps/qgis/./python/plugins\processing', 'C:/PROGRA~2/QGISCH~1/apps/qgis/./python', u'C:/Users/HEC/.qgis2/python', u'C:/Users/HEC/.qgis2/python/plugins', 'C:/PROGRA~2/QGISCH~1/apps/qgis/./python/plugins', 'C:\PROGRA~2\QGISCH~1\bin\python27.zip', 'C:\PROGRA~2\QGISCH~1\apps\Python27\DLLs', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\plat-win', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\lib-tk', 'C:\PROGRA~2\QGISCH~1\bin', 'C:\PROGRA~2\QGISCH~1\apps\Python27', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\site-packages', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\site-packages\GDAL-1.11.0-py2.7-win32.egg', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\site-packages\PIL', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win32.egg', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\site-packages\python_dateutil-2.2-py2.7.egg', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\site-packages\pytz-2014.2-py2.7.egg', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\site-packages\win32', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\site-packages\win32\lib', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\site-packages\Pythonwin', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win32.egg', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\site-packages\six-1.6.1-> py2.7.egg', 'C:\PROGRA~2\QGISCH~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', 'C:\Users\HEC\.qgis2\python\plugins\mmqgis/forms',     'C:\PROGRA~2\QGISCH~1\apps\qgis\python\plugins\fTools\tools']

>

The plugin for Cartogram Creator for qGIS is downloaded from here: https://github.com/cfarmer/cartogram-plugin/tree/master


texts: Are there any tools for editing Sql Server spatial types in Silverlight?

texts: I see from this post that I can tweak RIA Services to transmit Sql Server spatial types.  What I'd like to do is edit these spatial types.  Does anyone know of any tools that would run in Silverlight that would support this?


texts: How to prevent vector layer zooming if base layer does not allow farther zoom

texts:
I am using different types of base layer in OpenLayers and drawing vector layers on them. I am having trouble with Google satellite base layers.

I have noticed that satellite view has different zoom levels in different places. In some places base layer is not zooming but vector layer is zooming.

How can I determine if farther zoom level is available for a certain base layer?

In the images bellow a corner is marked with red. Notice the change position of the vector layer

Vector at zoom level 19

Same vector at zoom level 20


texts:
I'm trying to mask a raster with a ring polygon using the function raster::mask, but only recognizes the outer edges of the polygon, can anyone help me with this?
I only need the pixels within the yellow area, but it returns pixel inside the outer edge of the polygon.

library(raster)

raster<-raster("my_raster_path")
buffer<-shapefile("my_shapefile_path") #yellow buffer

output<-mask(raster, buffer)
writeRaster(output, "filename") #output grayscale raster


texts: Conversion between unknown map projections?

texts:
So I have two maps

One that divides regions into small provinces (represented with "ID"s)
Another showing which country controlled each part (shown with colours)

I want to overlay them to see which provinces would be controlled by each country.

The issue is that the maps have different projections. The one with the provinces is from Hearts of Iron IV and since I'm in Europe, it essentially is using the Miller Projection.

But I'm not sure which projection the other one is using. This was my attempt to overlay the two:

My question is: What is the best way to do this? Is there a program for it? Should I try even thinner strips? Maybe a mathematical approach is best?

Here's the second (non-province-containing) map if anyone thinks they'd be able to identify the projection.


texts: Viewport bounds online update. OpenLayers

texts:
I am trying to extract the bound box coordinates of viewport (the box that is displayed). I want to know the coordinates all the time and updated when zoom in/out or dragged. I cheched that 'getExtent' function returns the Bound coordinates, but I am not able to update the values. The code used is the one below (the script).

<script>
    var map;
    function Initialize() {
        LoadMap(40.36, -3.66, 6);
    }
    //***************************************
    function LoadMap(lat, lon, zoom) {
        var fromProjection = new OpenLayers.Projection("EPSG:4326");
        var toProjection = new OpenLayers.Projection("EPSG:900913");
        var lonLat = new OpenLayers.LonLat(lon, lat).transform(fromProjection, toProjection);

        var options = {
        projection: "EPSG:900913",
        center: lonLat,
        units: 'km'
        };

        map = new OpenLayers.Map("mapDiv", options);
        var osmLayer = new OpenLayers.Layer.OSM(null, null, zoom);

        map.addLayers([osmLayer]);
        map.addControls([
        new OpenLayers.Control.PanZoomBar(),
        new OpenLayers.Control.LayerSwitcher()
        ]);

        //Zoom the map
        map.zoomTo(zoom);

        $("#Bounds1").text("    ->     "+map.size);
        $("#Bounds2").text("    ->     "+map.getExtent());
        $("#Bounds3").text("    ->     "+map.getExtent().transform(toProjection,fromProjection));
        $("#Bounds4").text("    ->     "+map.calculateBounds());
    }
</script>

In this example, the map loads on Madrid and displays the coordinates of bound correctly although it seems to be wrong because  the 'zoomTo' function rapidly changes the view to afar view of all Spain.

So, how can the values on #Bounds labels can be updated? If it is not possible to do it in real-time, they can be updated each second or similar?


texts: Using PyQGIS, ArcPy, NumPy, SciPy Together in Windows Command Prompt

texts:
I ran across this slight issue when I was looking to access python installations from ArcGIS, QGIS, Anaconda (cited because SciPy & NumPy come with it), and version 2.7 vs 3.  Sometimes, I find it easier to just pop into Windows' Command Prompt to do a little bit of quick work instead of opening up a GIS program or python IDE in order to access the interactive prompt.
Instead of navigating to each folder, setting shortcuts, or switching out the PATH environment variable, I wanted a different way to access each python installation via the Command Prompt.
Edit
GIS professionals who use python for their work end up having to access each installation manually in order to access the features held within.  Instead of waiting on something like ArcMap to load, it is generally much quicker to hit your winkey-r + cmd and type python into your Command Prompt to start up the python console.  If you want to use ArcPy to reproject a shapefile, you need to open your python console, which is linked to ArcPy.  If you want to use a tool from PyQGIS because you don't have the appropriate ESRI license, you need to find the python.exe that is linked to those modules.  This can be a pain when you're trying to do work in ArcPy, PyQGIS, NumPy, or SciPy without having all the libraries installed cohesively as one python installation.
Typically, you're bound by your Path variable to which python installation is your default while in the Command Prompt.  This question and the answers following look to address how GIS professionals might be able to utilize multiple python installations together in a single environment.  Python IDEs target a single python installtion and require changing settings in order to utilize different ones.  For instance, my PyScripter install targets my ArcGIS installation of Python.  So, if I want to utilize my PyQGIS modules, I need to open up QGIS or change the settings in PyScripter.
When writing modules/full scripts, it makes sense to target one individually.  That way you might get some sort of code completion happening.  However, when you're just utilizing some tools in each installation to manipulate a data set, it would be much nicer if you can quickly flip python consoles/installations within a single screen aka the Command Prompt.
To those voting to close/off-topic...
Although this obviously applies to anything the python language is attached to, for the poor GIS user who can't afford to buy ESRI extensions for data analysis tools, it serves as a way for them to simplify their life.  As some of the answers have shown, being able to customize this process more, e.g. add arguments to .cmd files, could potentially be of great use if someone is running a custom process.  There is potential for a more savvy GIS/computer user to be able to set up custom .cmd files to quickly run processes without having to go through the fuss of dealing with different IDEs/GIS software installs etc..
The only downside to running a python console within the Command Prompt is the loss of code completion.  However, if you run something enough, such as arcpy.AddField_management(), the code completion just gets in the way.


texts: Aggregating taxi data to census block data based on number of taxi trips?

texts:
I am doing a study on taxi demand using block level census data and taxi drop off and pick data in the form of Latitude and Longitude around the united states.

I am hoping to create a map that shows block level counts for taxi demands based on pickup and drop data.

I am trying to aggregate taxi data at the block level then visualize it using arcmap or other visualization tools.

Due to the massive size of my data for the taxi data, I am trying to focus first on Delaware to test the waters.

I have two databases loaded up im MySQL with the Census data and the Taxi data.

How do I go about aggregating data at the census block level?

I have my script for the county level but I do not know how to do it for the block level.

Maybe we can modify my county level script to a block level script.

select * from [dbo].[BLOCK2010]
where STATEFP10 = 10

select geometry::UnionAggregate(Geom) as geom from [dbo].[BLOCK2010] where STATEFP10 = 10
Group By COUNTYFP10

declare @DE geometry = (select geometry::UnionAggregate(Geom) as geom from [dbo].[CB2010_PH] where STATEFP10 = 10)
--select @DE
SELECT [GEOM_BeginTrip] from [dbo].[TAXI_DATA] where [GEOM_BeginTrip].STIntersects(@DE) = 1 union all select @DE


texts: Seeking alternative site to GLC30/globeland30?

texts:
In 2014, China released a global land-cover data at 30m resolution, called GLC30 or globeland30.

I used it with success in some areas. But for some months now when I try to access the main download site www.globallandcover.com, I got a time-out error.
glc30.tianditu.com is still accessible, though.

Does anyone know an alternative download site?


texts: Converting a 2x2km point grid to a coincident 1x1km grid

texts:
I needed to convert a 2x2km point grid (specifically the LUCAS grid) into a coincident 1x1km grid.

Is there a "simple" way of doing this in QGIS?


texts: Looking for free/paid RSS to GeoRSS converters

texts:
I am looking for an RSS to GeoRSS converters. I was able to find only one converter i.e. GeoNames (it's crashing sometimes) and Meta-carta (looks like closed)

So does anyone knows free/paid RSS to GeoRSS converters?

Software: ArcGIS Suit 10.1 & ArcGIS JSAPI 3.3/3.4


texts: Open the IGxDialog in MS Access 2010 using VBA

texts:
I hope this is the right forum for this question.  I have a program in MS Access that allows users to input varius aspects related to land treatments including the spatial information (point, line, or polygon).  As part of the program the user can select a feature (shapefile) and the feature is automatically copied to the correct location, attributes attached, and loaded to a geodatabase.  Currently I use the standard file dialog Set fDialog = Application.FileDialog(msoFileDialogFilePicker) to accomplish this but as ESRI moves more towards geodatabases this is no longer an option as you can only select shapefiles (.shp) with the standard file dialog.  Selecting geodatabase (.gdb) does not produce the features within the geodatabase.  I want to switch over to ArcMaps file dialog (IGxDialog). The code I found (see below) opens the file dialog but then immediately crashes and shuts down Access before I make a selection. Am I missing something, do I need to adjust a setting or work with something like ArcObjects in Python, or why is it crashing my Access program?  Help I'm just an ecologist and in way over my head.

Private Sub Test_Click()
    'Set variables
    Dim pGxDialog As IGxDialog
    Dim pFilterCol As IGxObjectFilterCollection
    Dim pCurrentFilter As IGxObjectFilter
    Dim pEnumGx As IEnumGxObject
    Dim GetFile As String

    '  Select Case aFilter
    Set pCurrentFilter = New GxFilterFeatureClasses
    aTitle = "Select Shapefile or Feature"

    Set pGxDialog = New GxDialog
    Set pFilterCol = pGxDialog
    With pFilterCol
        .AddFilter pCurrentFilter, True
    End With
    With pGxDialog
        .Title = aTitle
        .ButtonCaption = "Select"
    End With

    If Not pGxDialog.DoModalOpen(0, pEnumGx) Then
        MsgBox "No selection : Exit", vbCritical
    End If

    GetFile = pEnumGx.Next.FullName
End Sub


texts: How to Modify and Combine rows in an attribute table?

texts:
I would like to modify this attribute table.

This is an example

First arrange column “ORG_FID” in ascending order
If the number of rows corresponding to the column “block Number” is same as well as “Lot Number”, then as per the column “ORIG_FID”, I want to add additional column as “ Combinetxt” and the combine text will arrange as per the “ORIG_FID” ascending order “N/F ONIS B. PUPO AND ALEIDA  M. PUPO, H/W” for both the rows.
Then I can dissolve the attribute table  on the basis on 3 columns. (“Block Number”, “Lot Number”, “Combinetxt”)


texts: ArcObjects interface to construct circles with segments instead of arcs

texts:
I am trying to see if it is feasible, using ArcObjects (VB.NET), to create a tool to construct circles comprised of line segments instead of arcs.  I've looked over the available interfaces here but do not see one for the task.

Is anyone aware of an interface to make this fairly straight forward?


texts:
I have a slight issue I just cannot solve. I use QGIS 3.6

I have multiple layers:

layer with multiple points with milage numbers;
a layer containing one or more polygons (let's call them "areapolygons").

I need to cut that (area) in pieces according to that milage points.

I tried the following, but that did not work:

create paths from the points layer
buffering the paths, so I get polygons (let's call them milage-polygons)
then clipping the milage-polygons with the "areapolygons"

I hope that'll give me areas of 100m. But I got some errors.

I need a sort of 'cutting' lines on the place of the milage points, to cut the polygon. How can I perform that?


texts: missing tables vertices_tmp and vertices_tmp_id_seq pgrouting workshop

texts:
I'm using OSgeo 6.5 dvd and studing pgrouting workshop for get shortest path.
I created database call "routing" and used "template_routing" for create tables but only 3 tables created. So then I executed following commands

# become user "user" (or run as user "user")
sudo su user

# create routing database
createdb routing
createlang plpgsql routing

# add PostGIS functions
psql -d routing -f /usr/share/postgresql/8.4/contrib/postgis-1.5/postgis.sql
psql -d routing -f /usr/share/postgresql/8.4/contrib/postgis-1.5/spatial_ref_sys.sql

# add pgRouting core functions
psql -d routing -f /usr/share/postlbs/routing_core.sql
psql -d routing -f /usr/share/postlbs/routing_core_wrappers.sql
psql -d routing -f /usr/share/postlbs/routing_topology.sql

then 7 tables created but few tables are missing when compare to workshop..

Workshop output;

List of relations
Schema |        Name         |   Type   |  Owner
--------+---------------------+----------+----------
public | classes             | table    | postgres
public | geometry_columns    | table    | postgres
public | nodes               | table    | postgres
public | spatial_ref_sys     | table    | postgres
public | types               | table    | postgres
public | vertices_tmp        | table    | postgres
public | vertices_tmp_id_seq | sequence | postgres
public | ways                | table    | postgres
(8 rows)

My output;

List of relations
Schema |       Name        | Type  | Owner
--------+-------------------+-------+-------
public | classes           | table | user
public | geography_columns | view  | user
public | geometry_columns  | view  | user
public | nodes             | table | user
public | spatial_ref_sys   | table | user
public | types             | table | user
public | ways              | table | user
(7 rows)

I got error like this

NOTICE:  CREATE TABLE will create implicit sequence "vertices_tmp_id_seq" for serial column "vertices_tmp.id"
CONTEXT:  SQL statement "CREATE TABLE vertices_tmp (id serial)"
PL/pgSQL function "assign_vertex_id" line 15 at EXECUTE statement
#########################
size of streets: 4093
size of splitted ways : 17328
finished

when executing this command

osm2pgrouting -file "data/sampledata.osm" \
                      -conf "/usr/share/osm2pgrouting/mapconfig.xml" \
                      -dbname routing \
                      -user user \
                      -clean

how do I get create missing tables


texts: Spatial Join in Python Field Mapping with Multiple Merge Rules

texts:
I'm performing a Spatial Join with a polygon (target) and a point (joinFeature)- which has 5 fields in it that I need to "sum". I have written a Spatial Join and (after many long days) gotten the field mapping to work! But this time is different because I am using a Merge Rule.

Here is my code, right now it breaks at the spatial join and outputs an empty feature class.

Also, this is in a Python Toolbox and uses Parameters to get the target and joinFeature. Like I said, I have gotten this to work before without a merge rule and with different fields...

        # Create a new fieldmappings and add the two input feature classes.
        fms = arcpy.FieldMappings()
        fms.addTable(target)
        fms.addTable(joinFeature)

        # Create the required FieldMap and FieldMappings objects
        fm_apt_poly_APT_ID = arcpy.FieldMap()
        fm_stud_ELEM = arcpy.FieldMap()
        fm_stud_INT = arcpy.FieldMap()
        fm_stud_MID = arcpy.FieldMap()
        fm_stud_HIGH = arcpy.FieldMap()
        fm_stud_TOT_STUD = arcpy.FieldMap()
        fms = arcpy.FieldMappings()

        # Get the field names for both original files
        apt_poly_APT_ID = "APT_ID"
        stud_ELEM = "ELEM"
        stud_INT = "INT"
        stud_MID = "MID"
        stud_HIGH = "HIGH"
        stud_TOT_STUD = "TOT_STUD"

        # Add fields to their corresponding FieldMap objects
        fm_apt_poly_APT_ID.addInputField(target, "APT_ID")
        fm_stud_ELEM.addInputField(joinFeature, "ELEM")
        fm_stud_INT.addInputField(joinFeature, "INT")
        fm_stud_MID.addInputField(joinFeature, "MID")
        fm_stud_HIGH.addInputField(joinFeature, "HIGH")
        fm_stud_TOT_STUD.addInputField(joinFeature, "TOT_STUD")

        # Set the output field properties for both FieldMap objects
        apt_poly_APT_ID = fm_apt_poly_APT_ID.outputField
        apt_poly_APT_ID.name = "APT_ID"
        fm_apt_poly_APT_ID.outputField = apt_poly_APT_ID

        stud_ELEM = fm_stud_ELEM.outputField
        stud_ELEM.name = "ELEM"
        stud_ELEM.aliasName = "ELEM"
        fm_stud_ELEM.outputField = stud_ELEM

        stud_INT = fm_stud_ELEM.outputField
        stud_INT.name = "INT"
        stud_INT.aliasName = "INT"
        fm_stud_INT.outputField = stud_INT

        stud_MID = fm_stud_MID.outputField
        stud_MID.name = "MID"
        stud_MID.aliasName = "MID"
        fm_stud_MID.outputField = stud_MID

        stud_HIGH = fm_stud_HIGH.outputField
        stud_HIGH.name = "HIGH"
        stud_HIGH.aliasName = "HIGH"
        fm_stud_HIGH.outputField = stud_HIGH

        stud_TOT_STUD = fm_stud_TOT_STUD.outputField
        stud_TOT_STUD.name = "TOT_STUD"
        stud_TOT_STUD.aliasName = "TOT_STUD"
        fm_stud_TOT_STUD.outputField = stud_TOT_STUD

        # Set the merge rule to sum and then replace the old field map in the mappings object with the updated one
        fm_stud_ELEM.mergeRule = "SUM"
        fm_stud_INT.mergeRule = "SUM"
        fm_stud_MID.mergeRule = "SUM"
        fm_stud_HIGH.mergeRule = "SUM"
        fm_stud_TOT_STUD.mergeRule = "SUM"
        fms.addFieldMap(fm_apt_poly_APT_ID)
        fms.addFieldMap(fm_stud_ELEM)
        fms.addFieldMap(fm_stud_INT)
        fms.addFieldMap(fm_stud_MID)
        fms.addFieldMap(fm_stud_HIGH)
        fms.addFieldMap(fm_stud_TOT_STUD)

        # Spatial Join Sample Code
        # SpatialJoin_analysis (target_features, join_features, out_feature_class, {join_operation}, {join_type}, {field_mapping}, {match_option}, {search_radius}, {distance_field_name})
        # Process: Spatial Join
        arcpy.SpatialJoin_analysis(target, joinFeature, output, "JOIN_ONE_TO_ONE", "KEEP_ALL", fms, "", "INTERSECT", "")
        arcpy.AddMessage("Just performed Spatial Join")


texts: Only first 10 layers of Geoserver WMS layers of enabled in ArcGIS

texts: I can see only first 10 wms layers of Geoserver is enabled while opening with ArcMap. Is there anyway, I can make all layers enabled?


texts: First report with a map has correct legend and scale. The second (other layer, new report) adds the legend of the first and has a bad scale

texts:
On QGIS 3.14 I've created a project.

I added a layer, created a map with quantile classification, then created a report for it. Added Legend, title, map, scale. Everything OK.

Then I did that :

I removed that layer.

I added another layer in order to produce a second report (I asked for a new report).

But even if only a single layer was displayed and checked in the Layer list :

Legend was displaying the legends of the current and the last map I did before.
Scale of that new map became false, on this report.

What could I have done wrong ?


texts: Give Universal BaseCommand a custom Icon

texts:
If I add ControlsMapZoomInTool like this:

AddItem("esriControls.ControlsMapZoomInTool", -1);

To a Basic Context Menu I get the magnifying icon in my context menu right next to the words "Zoom In".

Is there a way to add an icon to a custom Universal BaseCommand in ArcEngine?

ArcEngine 10, C#, VS2010


texts: Saving QGIS view extent in Shapefile or other exchange format?

texts:
I want to know how its possible to save the actual view of the map. It would be great to open it in QGIS or other GIS software the last saved view is showing when the file is opened. How can I load the saved shapefiles with the actual saved view?

When I open the last saved file, the actual view is gone, and I need to zoom again to the wanted map...


texts:
I have a shapefile with geometries of different regions. Also I have a list of coordinate points(field_lonlat), I want to consider only those points of the list which are outside the polygon geometries using python.
Below is my code to check it but trying to find a simpler version of it.
from osgeo import ogr
from shapely.geometry import Point
from shapely.prepared import prep
from shapely.wkb import loads
import json
g = ogr.Open("Polygons.shp")
points = []
for layer in g.GetLayer(0):
    field = loads(layer.GetGeometryRef().ExportToWkb())
    fieldP = prep(field)
    for x,y in field_lonlat: //field_lonlat is my list of points
        item = Point(x, y)
        if fieldP.contains(item):
            pass
        else:
            points.append((x,y))

But my points list checks for a single polygon once and hence considers only those points which are not in it but which might be inside an another polygon.


texts: Batch convert GeoTIFF's to Garmin Custom Map compatible KMZ files

texts:
I have about 600 GeoTIFF map files geo-referenced using New Zealand Transverse Mercator 2000 (NZTM2000). Each file is about 140 mb in size. I'd like to convert them to Garmin Custom Maps compatible KMZ files so people can download them and use them with Garmin GPS units.

Does anyone know a way to batch process these, preferably converting the tiff image to a much smaller jpeg image as part of the process? I'd prefer a tool that is free or donation based.

I've tried the method gdal_translate method provided as this answer https://gis.stackexchange.com/a/20936/1150
but unfortunately the KMZ files produced don't seem to be compatible with Garmin devices from the user feedback I've had. An example file produced using gdal_translate - https://dl.dropboxusercontent.com/u/48119/CB12-test.kmz

OkMap produces valid KMZ files but it lacks a useful batch processing option. OkMap also produces KMZ's that are much smaller that those produced through gdal_translate which is a nice bonus. Perhaps someone knows a technique / workaround I could use to automate OkMap not using the command line?


texts: Fields Read-only after Import from XML Workspace Document

texts:
I have a process that I follow where I import from XML into a Personal Geodatabase (PGDB), do some checks and updates, and then load into an enterprise GDB.

Now after I've imported from XML I currently have to copy/paste all the feature classes into another PGDB in order for the attribute fields to be editable.  If I try to edit them in the first PGDB, after they've been imported from XML, I just get messages along the lines of "field not editable"
When I right-click > Properties on one of the fields, the Field Properties dialog shows the field as being Read Only, and the option is disabled so I don't seem to be able to change it.

So, how can I either a) make the Import XML Workspace Document import as not read-only, or failing that b) make all my fields writable after doing the import, without having to copy all the feature classes first?
Running ArcGIS 10.7.1.  I have also tried doing the import into a File GDB with the same result.


texts: Different values for NDVI from different techniques

texts:
I'm trying to calculate NDVI from Landsat 5 imagery. I'm using ArcMap 10.2.2. I know that the formula for NDVI is (VIS-INFRA)/(VIS+INFRA), where in Landsat 5 VIS= band 3 and INFRA= band 4.

When I plug this into the raster calculator I get a raster with a max of 0.8, a min of -0.4 and a std of 0.1.

The NDVI tool in the Image Analysis window seems convenient (and would save a minute or two for each calculation if it worked properly), so I tried using that as well. I included a picture of how I'm using the image analysis window. I first designate the correct bands in the options menu, then highlight the proper bands and click on the NDVI button.

The raster I get from this tool has a max 0.98 a min of -0.98 and a std of 0.18. This histograms of both NDVI rasters look completely different.

Am I using the image analysis tool incorrectly? Or which approach in this case would give me the true NDVI values? I am leaning towards the NDVI values I calculated using the raster calculator as being correct, since an NDVI of 0.98 seems unreasonably high for my study area (urban).

Perplexingly, when I use the identify tool and pick any random location, I get the same NDVI value to six decimal places for both rasters. I'm not sure how that can be, yet to have the NDVI raster made through Image Analysis report a min, max, and std very different from the other raster.

Has anyone else run into a similar situation?


texts:
I am trying to calculate percentiles with gdal_calc.py using my command line.
I have a raster stack (11bands) and I would like to calculate the 0.85 percentile for each pixel over the entire stack (time).

Is there a way to do this?

Right now, I am using the numpy bindings of gdal_calc:

gdal_calc.py -A stack.vrt allBands=A --calc='nanpercentile(A.astype(int16),85,axis=0)' --outfile out.tif

I am trying to use the nanpercentile function from numpy. But I am not sure which parameters to use for axis and keepdims: https://docs.scipy.org/doc/numpy-dev/reference/generated/numpy.nanpercentile.html#numpy.nanpercentile

I am getting an output with this call, which however looks entirely wrong (the pixels look warped - it seems like the axis parameter is wrong).


texts: How to create dynamic interpolation maps using CartoDB?

texts: I can make interpolation maps in different time sets and save it as Rasters. How can I make them dynamic using CartoDB and its time slider?


texts: Can I pass multiple features to a Python Caller and use ArcPy in FME?

texts:
I've got a question about integrating ArcPy in FME. It would appear that FME won't let me to read in multiple features and process them using ArcPy functions through PythonCaller. This can be easily done outside FME with one line of code such as:

arcpy.LocateFeaturesAlongRoutes_lr(feats, rts, rid, tol, tbl, props)

What I would hope to achieve with PythonCaller is that it can be used just like one of the transformers like clipper and buffer. Not sure anybody else has experienced the same issue?


texts: Analyzing Raster Data From Vector Layer Buffer

texts:
I have satellite images of night-time lights all over earth. I also have a vector shapefile of country boundaries including the largest cities on earth. I want to use buffers to capture a certain amount of light pollution emanating from cities that is viewable from the raster image.

How do I use buffers to analyze pixels (or some aspect of the raster data) within the buffer area? I assume I will want to assign some value to a certain pixel and then total those values by using the buffer to analyze all data within the buffer.

I have attached a snip of what I am working with. I am very new to GIS still, so bear with me.

I am using QGIS 2.0.


texts: Getting Layer Properties of QGIS on standalone application using PyQGIS API

texts:
I am using PyQGIS API for developing the standalone application to load the certain VectorLayers. I used the QgsMapCanvas to load the map and QgsLayerTreeView to show the layer tree view using code below-

    # Layer Tree View
    self.layerTreeRoot = QgsProject.instance().layerTreeRoot()
    self.layerTreeModel = QgsLayerTreeModel(self.layerTreeRoot)
    self.layerTreeView = QgsLayerTreeView()
    self.layerTreeView.setModel(self.layerTreeModel)

    self.layoutWidgetLayout = QHBoxLayout()
    self.layoutWidgetLayout.addWidget(self.layerTreeView)
    self.layerWidget.setLayout(self.layoutWidgetLayout)

Now I want to style the layer and change some behaviors of the style using Layer Properties. How can I get that LayerPropertiesWindow by context menu of Layer as below.


texts:
I have visited, but do not remember where, Google maps for the Moon and Mars (I think). I am wondering if anyone knows what you have to do to generate the necessary GIS data (coordinate system, etc) for a completely new planet complete with aerials?

I would like to take other worlds and generate a Google Map complete with topography and draggability and more using the Google Maps API but I cannot find any information on how to do this online.

I would like to do this for both planets that are real and for fictional worlds. Think Narnia as a Google Map.


texts:
I've got a Layer out of an Shapefile with a lot different Attribute.
I colored them using properties->Style->categorized with random colours.
Now there are about 800 different colours. It looks something ike this. (It's in german).

So now if, let's say, I only want to see the second colour "abita Energie Ottenberg GmbH" on the map, the only way to do that is to uncheck all 800 other colours manually. I can't select more to uncheck them together.

Is there an easy way to do that? I hope you can help me.


texts: Two polylines to one singlepart polyline when endpoint share 3 vertices

texts:
I'm trying to merge/dissolve two polylines to create a singlepart polyline. The problem is that the two lines are a straight line with a loop (circle) at the end, meaning they share 3 vertices at the end point and so dissolve does not create a singlepart line. I need a line that goes towards the loop and counterclockwise around it ending where it meets the start of the loop and so far that only way I've found that works is to redigitize it manually.

The lines are road + cul de sacs.

Is there a tool or a way to make dissolve do this properly?


texts:
I followed this tutorial to set up Postgresql 9.1 and PostGIS 2 on my Ubuntu Server 12.04. I also installed pgAdmin III on my Windows 7 client. I set up a database landarch and user chadmwade, granting all privleges. I am able to connect to the postgres and the landarch databases with pgAdmin III and in QGIS.

I am getting the following error when attempting up add this shp file to the database. Is someone able to guide me in addressing this error? I haven't had much luck with Google.


texts: QGIS 2map plugin issue - no map tiles provided for QGIS 3x (3.8)

texts:
Following the example described here:

https://www.qgistutorials.com/en/docs/web_mapping_with_qgis2web.html

I am trying to do the same with QGIS 3.8.

The process in QGIS is simple, after seeking the following steps, being spread out in different places as the program has been update.
Example here:
Field editor in QGIS 3.0 no Text edition button

But the major problem for me is transfer this data into APIs provided (OpenLayers, Leaflet, Mapbox). After finish my process I get my data, but without any map tile...

what I can see already in the preview mode:

Thenceforth the map is blank even after export:

what I believe is caused by lack of option available in QGIS 2x, where user could select the map tilelayer at this step:

which diesn't appear in the QGIS 3.8.

Could anybody let me know, where can I find this option? I would like to have a full map exported not only the GeoJSON (layer) data.


texts: Create a layer, edit and save. Where are my features?

texts:
I created a layer (e.g. point) in QGIS with QgsVectorLayer from memory:

mem_layer = QgsVectorLayer("point?crs=epsg:32630&field=id:string","TempFile","memory")
QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

and save it with writeAsVectorFormat:

error = QgsVectorFileWriter.writeAsVectorFormat(mem_layer, "C:\\Temp\\TempFile", "CP1250", None, "ESRI Shapefile")

This all works well, but any features I add are not saved (only after save, remove, load again).

Is there anyone who can help me?

Here some more info. It is a tiny plugin to create a predefined SHP file:

(fileNameText,fileNameBase) = QInputDialog.getText(None,"QGIS in LUPMIS","Enter name of new point file:")
mem_layer = QgsVectorLayer("point?crs=epsg:32630&field=id:string&field=caption:string&field=style:integer",fileNameText,"memory")
QgsMapLayerRegistry.instance().addMapLayer(mem_layer)
temp1 = "C:\\LUPMIS\\LUPMISTemp\\" + fileNameText
error = QgsVectorFileWriter.writeAsVectorFormat(mem_layer, temp1, "CP1250", None, "ESRI Shapefile")
if error == QgsVectorFileWriter.NoError:
    iface.messageBar().pushMessage("LUPMIS","Point layer saved as file " + fileNameText + ".shp",QgsMessageBar.INFO)

From here, users should continue with manual entry through toolbar icons: Toggle Editing - Save Layer Edits - Add Feature, etc


texts: Using QGIS I am having issues in uploading a WFS, representing a geological map of Italy: when I upload it, the CRS sets automatically to EPSG:4326, ending up next to the Arabic peninsula. I used reproject layer, setting it to my project's CRS, EPSG: 32632. However, the new layer ends up in the same (wrong) spot. I tried to export and save as, then reproject layer, but it still ends up in the wrong spot. Is there something I am missing?


texts: I've come across code which specifies UTM=(32, 'U'). Is there a common convention around UTM I can understand this through, or is there not enough information here to determine the geographic location?


texts: Multiple Output for Zonal Statistics as Table

texts:
I am trying to merge NDVI data from 2005 to 2018 (two .tif per decade) with a region land use.

For this, I have used "iterate rasters" with 36 of .tif I'll be using, as the 'Input value raster' of the "Zonal Statistics as Table".
In the 'Input raster or feature zone data', I used a .shp with the regions I want.
Currently, it "works" manually and if the output table is fixed, but it overwrites the other generated tables (of course).
Using the %Name%_new however, I had errors (000354, 999999 and 010233) and nothing was generated.
I followed the steps in this link to solve my variable problem : gis.mtu.edu/wp-content/uploads/2013/04/Using-ModelBuilder-to-batch-process-files.pdf.

Do you have an idea of what to do? (in order not to overwrite again?)

The zonal statistics table content:

EDIT :
Bellow is my last modification:

I do not have errors but the out put is still overwritten. I feel like "iterate rasters" do not do his job properly in the model.


texts: Conversion of netcdf file to csv file using R programming

texts:
I have NCEP/NCAR reanalysis data in netcdf format. I want to convert it into csv file. I have a script to convert netcdf to csv. But since the data in netcdf file is too large which exceeds number columns in csv file. How can i split the data in nc file to multiple sheets of csv file.

Script that i have used to convert netcdf file to csv is as follows:

library(ncdf4)
getwd()
workdir <- "I:\\NCEP\\"
setwd(workdir)
ncin <- nc_open("X27.251.228.22.106.3.5.59.nc")
print("ncin")
dname <- "hgt"
lon <- ncvar_get(ncin, "lon")
nlon <- dim(lon)
head(lon)
lat <- ncvar_get(ncin, "lat", verbose = F)
nlat <- dim(lat)
head(lat)
print(c(nlon, nlat))
t <- ncvar_get(ncin, "time")
tunits <- ncatt_get(ncin, "time", "units")
nt <- dim(t)
tmp.array <- ncvar_get(ncin, dname)
dlname <- ncatt_get(ncin, dname, "long_name")
dunits <- ncatt_get(ncin, dname, "units")
fillvalue <- ncatt_get(ncin, dname, "_FillValue")
dim(tmp.array)
title <- ncatt_get(ncin, 0, "title")
institution <- ncatt_get(ncin, 0, "institution")
datasource <- ncatt_get(ncin, 0, "source")
references <- ncatt_get(ncin, 0, "references")
history <- ncatt_get(ncin, 0, "history")
Conventions <- ncatt_get(ncin, 0, "Conventions")
# split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth = as.integer(unlist(tdstr)[2])
tday = as.integer(unlist(tdstr)[3])
tyear = as.integer(unlist(tdstr)[1])
chron::chron(t, origin = c(tmonth, tday, tyear))
tmp.array[tmp.array == fillvalue$value] <- NA
length(na.omit(as.vector(tmp.array[, , 1])))
m <- 1
tmp.slice <- tmp.array[, , m]
lonlat <- expand.grid(lon, lat)
tmp.vec <- as.vector(tmp.slice)
length(tmp.vec)
tmp.df01 <- data.frame(cbind(lonlat, tmp.vec))
names(tmp.df01) <- c("lon", "lat", paste(dname, as.character(m), sep = "_"))
head(na.omit(tmp.df01), 20)
csvfile <- "cru_tmp_1.csv"
write.table(na.omit(tmp.df01), csvfile, row.names = FALSE, sep = ",")
tmp.vec.long <- as.vector(tmp.array)
length(tmp.vec.long)
tmp.mat <- matrix(tmp.vec.long, nrow = nlon * nlat, ncol = nt)
dim(tmp.mat)
head(na.omit(tmp.mat))
# create a dataframe
lonlat <- expand.grid(lon, lat)
tmp.df02 <- data.frame(cbind(lonlat, tmp.mat))
options(width = 110)
head(na.omit(tmp.df02, 20))
csvfile <- "cru_tmp_2.csv"
write.table(na.omit(tmp.df02), csvfile, row.names = FALSE, sep = ",")


texts: Drawing camera lines until frustum rectangle using Cesium?

texts:
I'm trying to draw a nice polygon like this to create my own sensor since we must to pay to have one from Cesium.

My goal is to make realtime viewshed analysis. For now I'm just making some tests. I don't know how to perform the viewshed "shadows". Will try to make it using the camera light and shadows.

To achieve this I need to 1) draw my own polygon like one of the picture but I confess I even know how to start. 2) use the Cesium Camera debug. This will draw a simple polygon simulating the FOV. Easy to do but...

The FOV polygon are going too far ( thru earth going to one galaxy far away...).

So how can I limit the "camera range" to the frustum (the rectangle near the camera, or the POV) ?

If someone have a more elegant solution I will gladly accept it.

My code until now...

var viewer = new Cesium.Viewer('cesiumContainer', {
    infoBox: false,
    selectionIndicator: false,
    shadows: true,
    terrainShadows: Cesium.ShadowMode.ENABLED,
    shouldAnimate: true,
    terrainProvider: Cesium.createWorldTerrain()
});

var scene = viewer.scene;
var ellipsoid = scene.globe.ellipsoid;

var instance = new Cesium.GeometryInstance({
  geometry : new Cesium.EllipseGeometry({
      center : Cesium.Cartesian3.fromDegrees(-100.0, 20.0),
      semiMinorAxis : 500000.0,
      semiMajorAxis : 1000000.0,
      rotation : Cesium.Math.PI_OVER_FOUR,
      vertexFormat : Cesium.VertexFormat.POSITION_AND_ST
  }),
  id : 'object returned when this instance is picked and to get/set per-instance attributes'
});
scene.primitives.add(new Cesium.Primitive({
  geometryInstances : instance,
  appearance : new Cesium.EllipsoidSurfaceAppearance({
    material : Cesium.Material.fromType('Checkerboard')
  })
}));

var west = -45;
var south = -24;
var east = -40;
var north = -20;
var homeLocation = Cesium.Rectangle.fromDegrees(west, south, east, north);

viewer.camera.flyTo({
    destination : homeLocation
});

var camera = new Cesium.Camera(scene);
//camera.DEFAULT_VIEW_RECTANGLE = homeLocation;
//camera.position = new Cesium.Cartesian3();
//camera.direction = Cesium.Cartesian3.negate(Cesium.Cartesian3.UNIT_Z, new Cesium.Cartesian3());
//camera.up = Cesium.Cartesian3.clone(Cesium.Cartesian3.UNIT_Y);
camera.frustum.fov = Cesium.Math.PI_OVER_THREE;
camera.frustum.near = 1.0;
camera.frustum.far = 2.0;

camera.setView({
    destination : Cesium.Cartesian3.fromDegrees(-45,-24, 15000.0)
});

scene.primitives.add(new Cesium.DebugCameraPrimitive({
  camera : camera,
  color : Cesium.Color.YELLOW
}));


texts:
I have a KML file (click link to display it in google maps) that works: https://drive.google.com/file/d/0Bx17tEcqrzavQUtzQzVqLTlkTVk/view

I have an embedded google map (works on my blog):

<iframe src="https://www.google.com/maps/embed?pb=!1m14!1m12!1m3!1d678888.4845987434!2d15.089529379165878!3d37.774381418009504!2m3!1f0!2f0!3f0!3m2!1i1024!2i768!4f13.1!5e0!3m2!1sen!2sus!4v1490795076461" width="600" height="450" frameborder="0" style="border:0" allowfullscreen></iframe>

How can I combine these two features so that the kml is displayed on the embedded map? I would like to avoid Java if possible.


texts: How to update GML layer without reloading local file?

texts:
I set up a GML layer with local GeoJSON file as data source:

veclayer = new OpenLayers.Layer.GML("vector", "localfile.json",
{
    format: OpenLayers.Format.GeoJSON,
    strategies: [new OpenLayers.Strategy.fixed()],
    styleMap: myStyleMap
});

Then later I change some of the attributes of GeoJSON data in the veclayer:

veclayer.features[i].attributes.someAttr = value;

After all features are updated, I want to redraw the layer. I see some examples that uses the following method:

function UpdateVecLayer(layer) {
        layer.loaded = false;
        layer.setVisibility(true);
        layer.refresh({force: true, params:{'key': Math.random()}});
    }

However, it doesn't work for me. I think it is because when I reload the layer, it just loads the local file again, and the attribute values are not changed. Can anyone tell me how to update the layer?


texts: QQIG/SAGA "merge shapes layers" algorithm removes underscores/pluses from layer names

texts:
I have a bunch of layers that are named layername_number+number that I want to merge. I need their names as source attribute so I used the "merge shapes layers" algorithm found under SAGA in the processing toolbox. While that does add the source layer name to the attributes of the merged layer (as "LAYER"), it removes the underscore and plus. Is there any way to prevent this?

The layers have no unique naming pattern (you could have name_0+100 or name_010+0, both resulting in name0100) so I can't reinsert the underscores and pluses afterwards either.


texts: Doing Math calculations within ArcMap Dynamic Text?

texts: Is it posssible to perform any math on ESRI ArcMap dynamic text?  For example I have the area of a SDE park layer's SHAPE.AREA, but it is in m2 because the layer is in NAD_1983_UTM_Zone_17N.  I want to show the area in hectares (m2 / 10000).  How do I do this?  Is it possible to add two fields together and display it in dynamic text?


texts: How to represent different colors based on a field with representation rule (ArcMap)?

texts:
I've got a single polygon feature class with hierarchy overlapping polygons (one-to-many - "tree"). The hierarchy is varying between 1-4 levels but can edited to be even more leveled in the future. The polygon feature class contains about 120 polygons with different NAME field values (though not completely unique).

I'm looking to distinct between the different polygons in two ways:

Buffer each polygon according to its hierarchy - the more basic it is (more a root in a tree) the more buffered.
Make a color ramp for the 'stroke layer' in the representation rule for each value in the NAME field.

I've achieved the buffering: I made  a double field and gave each feature a different value according to its hierarchy level; making an override to a buffer effect with that field.

My question is this:

How to generate different colors for the 'stroke layer' for that specific rule, according to the NAME field?

I've tried to follow this post (http://mappingcenter.esri.com/index.cfm?fa=ask.answers&q=2206) but the solution they suggested is to make a BLOB field and then manually change the colors for each feature with the representation toolbar, which can be quite a hassle.


texts: QGIS, Arabic names reversed (left to right instead of right to left)

texts:
I am working in a map with mixed names (Arabic and English). Data comes from OSM, a PBF export from BBBIKE, imported to Spatialite using OGR2OGR. I am adding the spatialite database using UTF-8 encoding.

Using the same display technique for both points and lines, for the POINTS layer I am getting the names displayed in Arabic as expected, without issues. However for the LINES layer, the names get displayed with the correct Arabic character but reversed (instead of being displayed for reading from right to left, are displayed in inverted order, as if they were meant to be read from left to right).

Here an example. Screen capture from the original OSM site of a section with a point feature (the mosque) and a line feature (the Brega Road, in red, at the lower right side of the image)

Then here is the result I get in GQIS:

Note that the mosque name gets displayed correctly in QGIS (font is slightly different, but the order is the correct one). On the other hand, the Brega Road below gets displayed inverted. That's it, should say "AQIRBLA QIRT" (for "TRIQ ALBRIQA", the name of the road in Arabic characters) but it says instead "TRIQ ALBRIQA" displaying the letters in the correct ORDER but in the wrong direction (as if it was meant to read like in English, from left to right).

Here below the capture of the content of the cells with the names. In both cases the content is the correct one, displayed correctly both in order and direction. The mosque name cell:

The road name cell:

Any idea on what could be wrong?


texts: How can I reproject a layer in a spatialite database

texts:
I'm writing an android mobile application and, for performance reasons, I want to reproject my layers to match the basemap's projection.

Is there a command to do that?