texts: Does SearchCursor reset work in ArcGIS 10.0?

texts:
In a script I have designed a function that is used multiple times on different search cursors to search through a whole table.

Function variables:

p - a list of points: actually a list of two-element lists with point coordinates, e.g. [[4710000.0, 3300000.0], [4715000.0, 3305000.0], [4710000.0, 3310000.0], [4705000.0, 3305000.0], [4715000.0, 3300000.0], [4715000.0, 3310000.0], [4705000.0, 3310000.0], [4705000.0, 3300000.0]]

c - a SearchCursor on a table with point coordinates

The function is supposed to search through the table with a cursor and return a list of points from the table that are also on the "p" list.

def ptsrch(p,c):
    list = []
    for i in c:
        currpt = [i.getValue("X_GRID"),i.getValue("Y_GRID")]
        if currpt in p:
            list.append(currpt)
    c.reset()
    return list

It works fine on the first use of each SearchCursor. Even though at the end of the function the cursor is supposed to reset, it doesn't.

In ArcGIS 10 Desktop help reset() is listed as a method: http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//000v00000068000000

I found a worrying post from 2008 on ESRI Forums:
http://forums.esri.com/Thread.asp?c=93&f=1729&t=250536

in which it is stated that "As far as I know, the Reset Method has never worked. I was told 2-3 years ago that it would be operational in a "future release". Apparently that time hasn't come yet!"

That one though refers to ArcGIS 9.3 or earlier.

I have also found a blog post: https://omniolocalia.wordpress.com/tag/arcpy-da-searchcursor/ which praises the possibility to reset cursors in ArcGIS 10.1, mentioning that it was previously impossible.

So, has anyone encountered a similar problem?

Instead of a SearchCursor I can of course load all tables to Python lists and keep them in memory, though I'm afraid of a memory overload - the total number of rows in these tables will be almost 400000 (right now I'm testing on a small chunk of data, so the total number of rows is 500).


texts: Correlation between raster values and vector point values

texts:
I will be using aerial photography/Object Based Image Analysis to separate plants in a field from their background, creating a polygon for each plant.

I will then convert these polygons to point on centroid and retain each polygon's area attribute, yielding DATASET 1.

I have a raster map of soil moisture of this field, DATASET 2, upon which I will overlay and georeference DATASET 1.

I will be normalizing the data so all raster pixel values in DATASET 2 and point values in DATASET 1 fall between 0 and 1.

I want to analyze the spatial correlation between the values of the points in DATASET 1 (plant area) collocated with values in DATASET 2 (soil moisture).

Can anyone suggest a technique to execute this? I have done a lot of research yet am finding it hard to find precedent to wrap my head around the problem.

I will be using ArcGIS 10.3 for Desktop to process the datasets.


texts: How to override shapefile geometry type in OGR

texts:
I'm new to ogr. I tried to copy a WFS layer to shapefile using CopyLayer but I get this error.

ERROR 6: Geometry type of `Multi Curve' not supported in shapefiles.
Type can be overridden with a layer creation option of
SHPT=POINT/ARC/POLYGON/MULTIPOINT/POINTZ/ARCZ/POLYGONZ/MULTIPOINTZ/MULTIPATCH.

How can I override the geometry type? Where do I apply the SHPT option?
from osgeo import ogr

url = 'https://data.gov.au/geoserver/brimbank-stormwater-pipes/wfs'
wfs_drv = ogr.GetDriverByName('WFS')
wfs_ds = wfs_drv.Open('WFS:' + url)
layer = wfs_ds.GetLayer()

driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open(shp, 0)
shp_layer = dataSource.GetLayer()
feature = shp_layer.GetNextFeature()
geometry = feature.geometry()
layer.SetSpatialFilter(geometry)
out_file = r'C:\stormwater-pipes.shp'
out_ds = driver.CreateDataSource(out_file)
out_ds.CopyLayer(layer, out_file)


texts:
I was istalled Mapserver:

sudo apt-get install libgd2-xpm-dev
sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
sudo apt-get update
sudo apt-get install cgi-mapserver mapserver-bin

and istall nginx:

sudo apt-get install nginx
sudo apt-get install spawn-fcgi

But what i gonna do now to Mapserver work in nginx? Can you give me a some links to tutorials?
OS: Ubuntu.

UPDATE

Okey. I found fastcgi config:

#!/bin/sh
MAPSERV_PORT=9001
USER=www-data
DAEMON=/usr/local/bin/mapserv
CHILDREN=5

# the -F switch of spawn-fcg does not work when the -n swich
# is set. using multiwatch instead
# see http://manpages.ubuntu.com/manpages/lucid/man1/spawn-fcgi.1.html

exec setuidgid $USER spawn-fcgi \
  -a 127.0.0.1 \
  -p $MAPSERV_PORT \
  -u $USER \
/usr/bin/multiwatch -f $CHILDREN $DAEMON

Now how to connect it to nginx?

UPDATE2

Now i add to the file /etc/nginx/sites-enabled/default follow lines:

   location /mapserver/ {
            include /etc/nginx/fastcgi_params;
            fastcgi_pass 127.0.0.1:9001;

            fastcgi_param SCRIPT_FILENAME /usr/local/run/spawnscript;
            fastcgi_param PATH_INFO /usr/local/run/ spawnscript;
   }

Where /usr/local/run/spawnscript path to my fastcgi config.

After i do:

/usr/lib/cgi-bin/mapserv -v

and get:

/usr/lib/cgi-bin/mapserv: /usr/local/lib/libgdal.so.1: no version information available (required by /usr/lib/libmapserver-6.2.1.so)
 MapServer version 6.2.1 OUTPUT=GIF OUTPUT=PNG OUTPUT=JPEG SUPPORTS=PROJ
 SUPPORTS=GD SUPPORTS=AGG SUPPORTS=FREETYPE SUPPORTS=CAIRO SUPPORTS=SVG_SYMBOLS
 SUPPORTS=ICONV SUPPORTS=FRIBIDI SUPPORTS=WMS_SERVER SUPPORTS=WMS_CLIENT SUPPORTS=WFS_SERVER
 SUPPORTS=WFS_CLIENT SUPPORTS=WCS_SERVER SUPPORTS=SOS_SERVER SUPPORTS=FASTCGI
 SUPPORTS=THREADS SUPPORTS=GEOS INPUT=JPEG INPUT=POSTGIS INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE

How i understand mapserver was installed correct.

Now i do:

sudo service nginx start

Go to localhos/mapserver/ but see only Welcome to nginx page. Whats wrong?


texts:
I am trying to build qgis without any change in source code which i have downloaded from github using clone and i have download OSGEO4W from qgis.org and the same is installed in c drive.

I have installed cygwin, cmake, OSGEO4w, mysys,python34, pyqtgpl, MinGW, QT5.5, qt creater, visual studio 2010 express.

Can any body share his cmake which has build successful.

I was able to configure cmake and generate. But I am getting 13 build failed in MSV 2010.
Im following the same path as give in source code procedure on qgis.org.

Can somebody help me out.


texts:
I have drawn a polygon in a shape-layer, but it's not perfectly circular. Is there a way to either

Draw a circular polygon with a given radius (in a shape-layer)
Make a polygon into a circle

I have tried to work around this by creating a point and styling it like my area, but then my point scales with the zoom levels, and it doesn’t look very nice (yes, I am aware of the different styles for different zoom-levels, but that's a lot of work compared to "just" making a circle-polygon instead)

I am using QGIS 3.2.3 with GRASS 7.4.1


texts:
I have a sh script that I want to use to process raster files.  eg

#!/bin/sh
# my first script,
# copyright, year, Author

r.mask -r

for file in A*chlor_a_4km ; do
  echo "map: $file"
    gdal_translate -a_srs "+init=epsg:4326" -a_nodata 65535 \
      -a_ullr -180 90 180 -90 -co "COMPRESS=PACKBITS" \
      $file ${file}_prep.tif
    r.in.gdal in=${file}_prep.tif out=$file --overwrite

done

I'm running from the OSGeo4W Shell, Windows 7, Grass 7.2.2, all freshly installed via OSGeo4W.

The script runs, but it produces eg 'r.in.gdal command not found' but then carries on to the next raster.

I expect this is something to do with the Path?  I've just installed this so can someone please send instructions on how to set the OSGeo4W Shell up for use.


texts: How do you specify the UTM zone in mapshaper?

texts:
I have a GeoJSON file containing boundary data for Norwegian administrative units called fylker. The original source of the data (GeoNorge) identifies the coordinate system as EUREF89 UTM zone 33, 2d.  I need it in WGS84.
I'm trying to use mapshaper to do the conversion.  I load up the file, open the console and execute the command -proj from=utm wgs84. It appears to work; the shapes displayed get slightly squashed, but I think that's normal.  After that I simplify the shapes a bit and export a fresh GeoJSON.
But when I open the new GeoJSON file in geojson.io, the shapes are offset quite a bit -- hovering in the ocean about half way between Norway and Iceland.  I'm not a GIS expert by any means, but it looks like it's interpreting the source data as if it were in UTM zone 30 or 31, instead of UTM zone 33.
I tried -proj from=utm33 wgs84, but it failed because it didn't recognize utm33 as a valid projection. The documentation doesn't mention this.
So how do I tell it that this data is in UTM zone 33?  Or is there something else I'm doing incorrectly?


texts: Integrating Google Earth historical imagery within Google Earth Engine

texts:
I'm planning on doing a high res. land cover classification of Southeast Asia for the last 20 years. I'd like to use Google Earth Engine Code Editor for the easy-access to Landsat data/ability to easily create time series/training data/etc., but also need to view high resolution historical imagery on occasion, e.g. the kind you can get via Google Earth Pro's historical imagery slider.

Is there some way to integrate the two within Earth Engine? Right now I'm just using coordinates to reference the same location between the two platforms, but  re-tying lat/lon again and again between the two platforms is inefficient and I need to streamline my training data collection process.

I know you can do Export.map.toCloudStorage(), but Google Earth web doesn't have the historical functionality....

Any suggestions on how to integrate the two?


texts: Removing duplicates from large LAS files

texts:
Are there any publicly and free (not necessarily open-source) applications that can remove duplicate points from LiDAR clouds?

I would like something on the lines of lasduplicate -unique_xyz, which seems to work fine, but LasTools are no longer free and have certain input limits above which the output is tainted by noise, missing information and it also seems to be incomplete.

Also, I've tried SAGA GIS's "remove duplicate points", but this option is almost impossible to work with for large sets (more than a couple of million and SAGA GIS hangs).

LATER EDIT: a Python script able to filter out such an input is provided in my answer to this question. You should force a garbage collect from time to time if you run into memory problems.


texts: How to use polygon to create raster landmask in QGIS?

texts:
I have a shapefile with polygons displaying the borders of the countries of Belgium, The Netherlands, and Luxemburg in EPSG:4326 (WGS 84).

How can I transform this into a raster landmask, showing values of 1 when inside the border and 0 outside the border? Additionally, the raster should have a 0.01° (~1.5km) resolution, with an extent between 49N-53.75N, 2E-7.5E.


texts: GeoPandas sjoin with extensive quantities

texts:
I have one geodataframe that contains polygons and corresponding land use category (e.g. agricultural, commercial, retail, residential, ...) and another that contains census tract polygons and corresponding population.  I want to bring the population data from the 2nd dataframe into the first dataframe.  Here is the code I'm using:

gdf_chi_pop = gpd.sjoin(gdf_chi, gdf_census, how='inner',op='intersects')

The problem is that gdf_chi_pop has a bunch of extra polygons that don't exist in gdf_chi, which is not what I'm wanting.  I just want a new dataframe with populations for all the polygons that already exist in gdf_chi.  What I also need to know is how geopandas differentiates between intensive properties (e.g. temperature, sunny days per year, population density, fraction of dogs that are poodles) and extensive properties (e.g. population, number of poodles, number of stores that sell spatulas) in polygons when doing sjoin.


texts: Seeking hardware recommendations for Mapnik Server?

texts:
I'm looking for some suggestions for a Mapnik-TileCache server.

How much hardware is really necessary?

The dataset I use is about 380mb in osm raw format. My actual database is a bit more than 1GB.

What would you suggest?

QuadCore
8GB Ram
SSD Harddrive for DB
Bunch of disks for Cache

What do you suggest?

The server will have frequent trafic and let's say 300 - 500 people will work on that mapset.


texts:
I am confused, whether  removing pseudo-nodes and line-simplification methods should be considered as the same concept.

Line simplification is one normally done using Douglas-Peucker algorithm such as with v.generalize tool in GRASS.


texts:
Where can I change the default labeling font "MS Shell Dlg 2" to a different one?

I don't want to change the default layout font (Options --> Generel Settings)

Maybe there is a solution to use the global variable, but is there an option to always activate the global variable for fonts?

The aim is to change the font here:

What i've already tried:

to find and change the font variable for style labelling

to change the default font from windows regedit settings

Also i found the option to set a default font in Print Composer, that works.

I use QGIS 2.18.17 LTR and QGIS 3.2.1


texts: How to improve LiDAR classification in GlobalMapper

texts:
I have a LiDAR point cloud that has a density of 16.4 points per m^2 and spacing of 0.2469m, and am trying to extract trees from it but the classification came with a lot of messy points in classes they shouldn't be in, so I am trying to reclassify it from scratch using GlobalMapper V21.

My problem is that when I try and classify the buildings and vegetation I am getting building roof outlines in the high vegetation, though the main building roofs do show up in the buildings classification.  I have tried many permutations of the classification filtering settings and have improved things a bit, but I still can not get these building outlines out of my vegetation layer.

I have read the GM help page on what the settings do but am still scratching my head as to the best values to use.  The ground classification help page had good descriptions and example pictures of results using differing values, but the building/veg classification pages don't have the same degree of explanation so i am finding it hard to trial & error my way to a good result.

What suggestions do you have for improving this so I get buildings in the building class and only trees in the high & medium vegetation classes?

Attached is a picture of the best result I can get showing the building outlines, but also the settings seem to miss some of the trees.


texts:
Using GDAL tools such as gdal2tiles.py, allows one to generate map tiles from a single image file.

However, the tiles generated by gdal2tiles doesn't seem to be like one used for google.maps.ImageMapType Overlay

How to generate tiles using gdal2tiles that can be used for google.maps.ImageMapType?


texts: Combining line segments into single line in Google Earth Pro

texts:
In Google Earth Pro I have a KMZ for a single route, but its segmented into many line segments.

I need to be able to export this KMZ as a single feature - one line for the whole route.

Am I able to combine all the segments into a single line in Google Earth Pro and, if so, then how?


texts: Viewing print text while running a script in ArcMap

texts:
Why can't I see print text when running a script in ArcMap, but can view it when run in the python window of ArcMap?

I have a number of scripts that have a print function but don't show when run from the Toolbox.


texts: Removing measure label from OpenLayers 3 measure?

texts:
I have implemented the measure example available on the OpenLayers 3 site (http://openlayers.org/en/v3.11.2/examples/measure.html) and I am currently trying to include a clear button that will clear both the drawn area / line and label from the map.

I've found the getSource().clear() option will clear the drawn shape from the map, but I am struggling to work out how I can clear the measurement label?

I have tried the map.removeOverlay() method with 'measureTooltip' but this doesn't work.

How can I also remove this label from the map?


texts: Listing all tables with no geometry using Oracle DB and ArcMap

texts:
I imported all my data (including attribute tables without geometry) to Oracle Database (default db : SYSTEM)

I connected it to QGIS and ArcMap.

I need tables that starts with RDF*. it shows in QGIS but not in ArcMap

The picture below is screetshot from ArcMap's catalog and as you see it shows only feature classes and random tables that i do not know where they are from

This is from QGIS Browser Panel and it has tables that i need, although it also has tables that I don't know where they come from.

I want my ArcMap catalog view as in QGIS's browser panel.
can anyone give me a hint?


texts: Creating barriers based on distance along geometric network

texts:
I have a task at work that requires that I select 1 mile of line (pipe) downstream from a specific point (manhole) in a network.  Each line in the network has a different length, so I need the sum of the length of lines to be around 5280 ft.  I have roughly 300 points within the network, so I'm looking for an automated way to do this.

So far, I've used the geometric network trace tool to trace downstream from each point (flag).  I need to add a barrier along the each downstream route at around 5280 feet.  Does anyone have an idea for how I can accomplish this?

I'm currently using ArcGIS 103.1 Advanced without a Network Analysis license.  I am also using the geometric network toolset to create my network.


texts:
I'm working on a project that requires custom styles for lines, similar to that in Mapbox's line-pattern-file example, but I keep running into an issue where my custom line symbology looks suboptimal on high resolution displays. Ideally, I would use SVGs instead of PNGs for this, but it seems that only rasters are supported by Mapnik for styling lines.

Although my use case deals with lines, more generally, are there any best practices or tips for working with raster fills and patterns on high resolution displays?

The documentation acknowledges this issue, but only briefly mentions that it's a good idea to create two versions of fills and patterns for this purpose. From my testing, creating a 2x version of pattern simply changes the scaling of the pattern and does not actually look higher resolution.


texts: Import of NAS-data, (norm-based data exchange interface)

texts: Is there any possibility to import nas-data in QGIS. NAS (norm-based data exchange interface) is a XML-based format. German cadastral-data is only available in this format.


texts: Are ArcMap Addins compatible with ArcGIS Pro?

texts: In other words, can ArcGIS Pro load an ArcMap AddIn?


texts: Replace Navteq maps by a free one in a car built-in GPS

texts:
My car (Opel Zafira) have a built-in GPS with Navteq maps (SD 900 NAVI).

I wonder if it would be possible to update or upgrade the maps to an open source / free maps (like OSM).

Or is it possible to convert the OSM maps into the Navteq format?

Anybody tried this?


texts: How to convert a microstation cell library into a format usable in ArcGIS?

texts:
I am in the process of standardising symbology across CAD and GIS and have decided to use the symbols we currently used for points in microstation. These are stored in a cell library.

Does anyone know how I can create a style or a font out of these cells and therefore use microstation symbols in ArcGIS?


texts: Help with symbol property editor in arcmap

texts:
I want to create symbols like below With different colouring in a polygon.

How do I do it?

any videos that can help me?


texts: How to count point features in polygon using ArcGIS API for JavaScript?

texts:
I have published in my ArcGIS Server two layers: A polygon one and a point feature one. Points are contained in polygons.

I would like to know how to count the number of points belonging/contained to a polygon using ArcGIS Javascript API. So, the functionality would be, clicking in one of the polygons and then, present somewhere an integer (in an Infowindow or textbox in sidebar) saying "There are X points in this polygon".

Could you provide any ideas or examples?


texts: Using .DXF Arc or Circle accurately in ArcGIS for Desktop?

texts:
I need to use a .DXF in either QGIS or ArcGIS.

However, if an Arc occurs in a .DXF drawing and if this file is opened in QGIS, I noticed that the geometry of the Arc is simplified. I did a test:

draw circles in LibreCAD and save as .DXF (I uploaded it here)
open in QGIS
convert to shapefile
edit shape

The circles, regardless of their size are approximated by 90 segments, so a 4-degree long circle segment is represented by a straight line. In my application this causes significant deviations (>25cm) in a canal design project...

I tried the same in ArcGIS. There, the simplification depends on the circle size - small circles have fewer segments, larger circles have more. So, in my case (for the kilometer long radius in canal bends) ArcGIS seems like the better choice.

As a side note: In ArcGIS the Arc of the .DXF is not simplified when visualized - the display of the .DXF seems correct when overlying with the exported shapefile. Within ArcGIS this would make at least snapping and digitization in higher accuracy possible.

I am a beginner in using .DXF files in GIS, but I feel this should be a popular problem. Though, after searching for a while I could not find a hint anywhere. I still fear I may miss something obvious... In case not, does anybody have a workaround for this?

I am using QGIS 2.10 (GDAL 1.11.2) and ArcGIS 10.2


texts: Are Planet Asset Activations Bundled or Synchronous?

texts:
I was wondering if planet asset activations are linked/bundled or synchronous for item types so for example
If I use --activate PSOrthoTile analytic, does it automatically activate analytic_xml and also analytic_dn? or would I have to run --activate PSOrthoTile analytic_xml and --activate PSOrthoTile analytic_dn separately?


texts: Where are ComReleaser & WebObject ManageLiftime in ArcGIS 10.2.1

texts:
I used to use ComReleaser in desktop, and WebObject in ArcGIS Server apps to release com objects - e.g. cursors and the like.

I am undertaking an upgrade from 9.3.1 to 10.2.1 of an app - and for the life of me I cannot find the ComReleaser and WebObject to help me manage com object lifetimes!

Are they just not needed anymore in 10.2.1?


texts: Estimating processing time for rasterize or polygonize tools in QGIS?

texts:
I often rasterize or polygonize large layers (flood modelling results).

Is there a way I can estimate how long the tool will take?

Sometimes, I'm not sure whether it will take 10 minutes or 10 hours or 10 days.


texts: How to display geotagged photos in leaflet

texts:
There is awsome plugin of displaying geotagged photos in leaflet
Photos on leaflet

I used the picasa photos implimentation
when I implemented the plugin for  my picasa public album photos link in this format
https://picasaweb.google.com/data/feed/api/user/MY-USER-ID/albumid/MY-ALBUM-ID?alt=json-in-script

my example failed even my link work perfectly and broadcast Json data

Can you show me how to use that plug in?


texts: Deleting traces after dissolving several polygons in ArcGIS Desktop?

texts:
I'm using ArcGIS 10.4.

I have a shapefile that has many catchments in it. I wanted to dissolve all these catchments to get just the boundary of the study area.

I used dissolve to do so

The problem is that I have found traces after I dissolved all polygons as below

How do I get rid of these traces?

DATA

You can download the shapefile from this link.


texts:
I want to intersect a Z-line with a non-Z-line to obtain the intersection point with the z-value from the Z-line.
However, it doesn't work for these two lines:
my_line <- st_linestring(matrix(c(403606.39, 403601.85, 403601.09,
                                  5306798.72, 5306802.81, 5306803.73),
                                nrow = 3),
                         dim = "XY")

my_zline <- st_linestring(matrix(c(403601.81, 403601.85, 403602.26,
                                   5306802.76, 5306802.81, 5306803.48,
                                   241.61, 241.73, 242.98),
                                 nrow = 3),
                          dim = "XYZ")

st_intersection(my_line, my_zline)

Result, without the desired Z-coordinate:
POINT (403601.85 5306802.81)

What am I doing wrong?


texts: Arcpy.da.SearchCursor to retrieve the min and max measure values from Polyline-M shape

texts: Is there a way to use a search cursor in ArcGIS to read the minimum measure and maximum measure values of a linear referenced route layer using arcpy? My goal is to build a route event using the start and end points of line segments but cannot figure out how to retrieve this data other than using the Identify Route Location Tool in ArcMap.  I would like to get that type of information but using python.


texts: Unable to store relative pathnames to data sources in ArcGIS Publisher

texts:
The Publisher extension in ArcGIS Desktop 10.5.1 isn't saving the relative data source for our enterprise geodatabase layers. The data sources are referencing the enterprise geodatabase directly rather than in the specified location set in the Data Packaging.

ESRI tutorial Link:
How To: Create a PMF that can be opened in ArcMap

ArcGIS Desktop 10.5.1 - Map Document Properties


texts: Federate GeoEvent Server or Register Portal to GeoEvent Server?

texts:
In GeoEvent Server, I have the option to register Portal. What is the difference if I federate GeoEvent Server to Portal and register Portal to GeoEvent Server? Or I should do both?

I am using ArcGIS Enterprise 10.5.


texts: How to perform If Statement in Field Calculator of ArcGIS for Desktop?

texts:
I have a shapefile containing two numeric fields ("Dist_1" and "Dist_2").  I want a Field Calculation that will populate an additional field ("Result") with one of three answers:

First if Dist_1 is greater than Dist_2 populate Result with 'True'
Second if Dist_1 is less than Dist_2 populate with 'False'
and Third if Dist_1 and Dist_2 are equal populate with "Equal"


texts: How does GPS receiver know estimated error?

texts:
I would like to plot a circle of possible error on the map in a GPS app.

From http://en.wikipedia.org/wiki/Dilution_of_precision_(GPS), we can see that GDOP is equal to "delta of output location" over "delta of measured data". How can we know the "delta of measured data" in order to calculate "delta of output location"?

If it is not possible, how do the receivers plot the circle of possible error on the map?


texts: How show some fields only when using WMSGetFeatureInfo?

texts:
How to show some important field only when using WMSGetFeatureInfo.I use this example.
When I click on maps it show all the fields of attribute table.What part should I modify in above mentioned code.
Pls help me?
Any idea wellcome


texts: ArcGIS Solar Radiation Tool - How is Local Solar Time/Local Time handled as input to the "within-day" time configuration?

texts:
I am scripting hourly interval runs of the Area Solar Radiation Tool, and need a hand understanding how to account for Local Solar Time (LST) versus Local Time (LT) within the tool.

First off, I understand the difference between Solar versus Local Time, and the practical/fundamental basis of the tool developed by Drs. Fu & Rich. For those who need some context for this question, please visit the following excellent resources:

  ESRI Help on Solar Radiation Tool / PVeducation.org's Solar
  Time / How Solar Radiation is Calculated

Here's the question: Considering the following paragraph from the first link:

  For within-day time configurations, the start and end times are
  displayed as solar time (units: decimal hours). Use the time
  conversion dialog box window to convert the local standard time and
  local solar time (HMS). When converting local standard time to solar
  time, the program accounts for equation of time.

If running this as a script, how would one construct their "TimeWithinDay" object with the above answer in mind?

I have looked into the TimeWithinDay help page, and here is the code provided as an example of how I might want to run the script:

# Name: TimeWithinDay_Ex_02.py
# Description: Execute AreaSolarRadiation using the TimeWithinDay object
# Requirements: Spatial Analyst Extension

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = "C:/sapyexamples/data"

# Set local variables
inRaster = "solar_dem"

# Create TimeWithinDay Object
day = 100
startTime = 0
endTime = 24
myTimeWithinDay = TimeWithinDay(day, startTime, endTime)

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute AreaSolarRadiation
outAreaSolar = AreaSolarRadiation(inRaster, "", 200, myTimeWithinDay, 14, 0.5,
                                  "NOINTERVAL", 1, "FROM_DEM", 32, 8, 8,
                                  "UNIFORM_SKY", 0.3, 0.5)

# Save the output
outAreaSolar.save("C:/sapyexamples/output/areasolartwd2")

Anyway, back to the question...does the above code consider the Start and End times "fed" into the tool, via the "myTimeWithinDay" object, automatically calculated into LST, or does it assume LT inputs and calculate accordingly? It notes accounting for the "Equation of Time", so what if I pass it LST times?

Given an answer to the above, is there a way to run the "dialog box window" converter for converting between LST and LT prior to (or during) the passing of Start and End times to the tool itself?

If not, I can calculate LST for LT intervals desired and feed LST times to the script to then have an hourly irradiance interval that has been calculated for LT, correct?

To help understand why I am interested in knowing how this works consider that I am wanting to base the transmissivity and diffuse fraction of empirical hourly values that were measured in Local Time. I need these two time series to match over a given hour of running the tool.


texts: PyQGIS: How to apply a color ramp on the symbol's outline when using QgsCategorizedSymbolRenderer?

texts:
I have polygon layer containing overlapping features (service areas for several years) and I want their symbology to

use a color ramp in QgsCategorizedSymbolRenderer applied for each year
use a transparent fill and apply their former fill color to the outline.

For the first part I successfully use this code:

categories=[]
for u in layer.uniqueValues(fieldIndex):
    symbol = QgsSymbol.defaultSymbol(layer.geometryType())
    category = QgsRendererCategory(u, symbol, str(u))
    categories.append(category)

renderer = QgsCategorizedSymbolRenderer("year",categories)
style = QgsStyle().defaultStyle()
ramp = style.colorRamp('Greens')
renderer.updateColorRamp(ramp)
layer.setRenderer(renderer)
layer.triggerRepaint()

How do I get each symbol's fill color after having applied the color ramp in order to make it the outline color for each symbol according to the ramp's gradient? Could there be a way to do so in the loop beforehand and then apply the ramp directly on the outlines?

Sidenote: I am really puzzled that using a feature's outline for symbology apparently requires some tinkering everytime...


texts: Drawing a Straight Line Between Two Polygons on a Map and Listing Which Polygons It Passes Through

texts:
I have a shapefile of United States counties.  I have two counties that I would like to highlight.  I would like to draw a straight line distance between them.  By straight line, I mean, if you were to lay down, say, railroad tracks that go straight.  I want to know which counties this line would pass through.

For example, suppose that I would like to draw a line between San Diego County, California, and Cook County (containing Chicago), Illinois.  I would want to know which counties this line would pass through.  For example, I imagine it would include counties in California, possibly Nevada, and so forth until Illinois, and obviously not include counties in states like Texas, Oregon, Florida, or Vermont.

I use ArcGIS and ArcMap 10.1.  I would very much appreciate any help that you could provide!


texts: What to do to keep Python code permanent in QGIS 2.14.2?

texts:
I have this code provided by @Joseph but lose it whenever I leave the QGIS and I have to load it on the Python Console whenever reopen QGIS.

What to do to keep this code in a permanent way?


texts: Why the execution of a query is very slow using PostGis?

texts:
I'm applying this Query to 8257800 point but it's taking 15 hours and it didn't return a results.

SELECT a.gid , b.gid, st_distance(a.geom,b.geom)
FROM shp1 a, shp2 b
where st_intersects(st_difference(st_buffer(a.geom,2000),st_buffer(a.geom,500)),b.geom) and abs(a.value - b.value)>400

Is there any suggestions to render this query fast?


texts: ESRI Javascript FeatureLayer Select by Attributes

texts:
I am trying to select a group of features in an ESRI javascript API featureLayer by attribute. The call errors out. My code is as follows. Here is the layer definition:

    var fieldDefs = [{
        "name": "county",
        "type": "esriFieldTypeString",
        "alias": "county"
    },
    {
        "name": "pop",
        "type": "esriFieldTypeInteger",
        "alias": "pop"
    },

        "name": "ObjectID",
        "alias": "ObjectID",
        "type": "esriFieldTypeOID"
    }
    ];

    var featureCollection = {
              "layerDefinition": null,
              "featureSet": {
                "features": [],
                "geometryType": "esriGeometryPolygon";
              }
            };
    var featureCollection.layerDefinition = {
              "geometryType": "esriGeometryPolygon";,
              "objectIdField": "ObjectID",
              "drawingInfo": {
                "renderer": {
                  "type": "simple",
                  "symbol": symbol
                }
              },
              "fields": fieldDefs
            };
     var featureLayer = new FeatureLayer(featureCollection, {
                  id: idString,
                  mode: FeatureLayer.MODE_SNAPSHOT,
                  fields: fieldDefs
                });

I then populate the layer using applyEdits and later query the layer using this function:

queryCounty = function(county_name) {
            var query = new Query();
            query.where = "county = '" + county_name + "'";
            //projectsOverlayMap.clearSelection();
            console.log("WHERE: " + "county = '" + county_name + "'");
            projectsOverlayMap.selectFeatures(query, FeatureLayer.SELECTION_NEW, function(results){
                alert(results.length);
            });

And calling it like this:

queryCounty("Jefferson");

When I do that I get the following error:

  "Error: FeatureLayer::selectFeatures - query contains one or more unsupported parameters

There are some decent on line examples of this function using geometry queries but none that I have found that use and attribute query. What am I missing?


texts: Calculating a point Lat/Lon from another Lat/Lon point angle

texts:
I need to make an algorithm:

:

I have a point O center of a circle of radius r, on this circle there is two point P1 and P2.

I know latitude and longitude of point O in degrees
I know the radius r in meters
I know the latitude and longitude of P1 in degrees

I want to know the latitude and longitude of point P2 on the same circle, I know the angle a in degrees between the line OP1 and OP2.

How to solve this problem?

The gps points O,P1 and P2 are situated in the south of France.


texts: Deleting overlapping inside same polygon using QGIS?

texts:
I have same polygon partially overlapped:

How to delete overlapping and merge it to same area?

I am using QGIS 2.4.0.


texts: Avoiding or dealing with "unable to save" errors in ArcGIS ModelBuilder?

texts:
Do you have references to best practice, or examples where you have recovered Models or can you share your expertise on how to avoid the "Unable to save model" error within ModelBuilder using ArcGIS 10.1 for Desktop?

Alternatively can you advise on situations to avoid which cause this error to occur.

There is very little help information on this topic.

There is a linked question here Do you use a proven workflow to edit and document ArcGIS ModelBuilder models?
 - that shares advice on how to track the building of models, so you always have a clean copy of the model to go back to. Comment there if you have comments about the process of tracking edits / changes to Models.


texts: Breaking long line segments into shorter ones in R using sf

texts:
I have a shapefile of river segments that I've imported into R. The segments range in length from 5-115km. I want to split the segments that are over 10km into 10km chunks. So I use st_segmentize. But this doesn't actually split a long line into multiple new shorter lines, which is what I want... Is there a way to do this? I think st_segmentize just adds vertices every 10km. So I just need to split the lines at the vertices.

Code so far:

library("sf")

rivers = st_read("./YukonMain4_copy.shp", layer = 'YukonMain4_copy')

seg = st_segmentize(rivers, units::set_units(10,km))

Re. my above comment

The issue is, my resulting lines have lengths mostly from 200-500m, and then there are random onesthat are 700000m long. Plus, this data set is of multiple rivers, so now I have lines connecting places that should not be connected. Any ideas?

Here is my code:

Import file, segmentize the data, and cast it into the point cloud
Create an index to combine segments by (hopefully so that they are approx 10km in length
Combine things back into lines.

rivers=st_read("./YukonMain4_copy.shp", layer='YukonMain4_copy')
seg = st_segmentize(rivers, units::set_units(10000,m))
seg2= st_cast(seg, 'POINT')

length = length(seg2$OBJECTID)/2
index = rep(1:length, each=2)
rivers.index=cbind(seg2, index)

to_line <- function(points) st_cast(st_combine(points), "LINESTRING") %>% .[[1]]
rivers.nest = rivers.index %>% group_by(index) %>% nest
segments = rivers.nest %>% pull(data) %>% map(to_line)%>%st_sfc(crs = 4326)
lines = rivers.nest %>% select(data)%>% st_sf(geometry=segments)


texts: Remove a custom contents view from ArcMap TOC in VB .NET?

texts:
I am using the ESRI sample, TOCLayerFilterVB2010, which implements a custom TOC view for map layers, as a prototype.  For my purposes, I need to be able to turn it on/off, meaning add and remove that view from the TOC horizontal icon bar programmatically.

Can someone explain how the custom contents view in this sample gets instantiated?  I find that if the dll is registered, it appears. The only way to get rid of it in any ArcMap document is to unregister it. Deleting the Normal.mxt template does not have any effect.

Ideally, I would like to have a command toggle that instantiates and destroys the view. I cannot figure out how to control its creation.


texts: GDAL error in QGIS application displaying VRT files

texts:
I am working with a standalone app build with QGIS 2.8.1 on a CentOs7 machine. I got a new DTED data set which raises errors "somewhere" in the app. The app doesn't read the .dt? files directly, but as a VRT.

Even in debug mode this is all I get on the console:

Warning 1: Longitude count (195) of column 193 doesn't match expected value
Process finished with exit code 139.

Just for info:

About QGIS states: compiled with GDAL/OGR 1.11.1, running with GDAL/OGR 1.11.2; Python is running with gdal 1.11.2
other DTED files work perfectly, readed directly or as VRT
reading a single .dt? file works
both, VRT and individual .dt? files work in QGIS application
The columns may be different in each run, and there may be different counts of warnings

Unfortunately yesterday evening CentOs was proud to tell, that QGIS has successfully updated to the current version 2.8.3. But using this version, and after downgrading back to 2.8.1, there is no difference. I cannot state, if this error would have occur before updating, because I hadn't used the same DTED data before.

What is the reason for this GDAL error, which instructions will raise this exception, and why is the behaviour of the different components (VRT/direct reading, QGIS/QGIS components) so different?


texts: Plugins installation failed because "network access is disabled"

texts: Any idea how to solve the problem with the plugin installation failed in qgis 3.0? I cannot run the installation for some plugins and it is because "network access is disabled".


texts: Why is a stored procedures much slower than the single query?

texts:
I wrote a normal query to snap the next way point from a pair of coordinates in EPSG4326, nothing complex. To optimize the performance a bit I made a boundary condition to minimize the points needed to be checked...

Works smooth so far - but then I wrote a stored procedure to make later a neat intergeneration... (quasi wrapper). And here is where my headache starts:

Single SQL: 650ms
Stored Procedure: 45xxxms
Single SQL is 70x faster than my storage procedure!

Code below...

CREATE OR REPLACE FUNCTION next_way_point(p geometry) RETURNS next_way_point_result AS $$
SELECT gid, source, target, the_geom, ST_distance(the_geom, ST_SetSRID($1,4326)) as dist
FROM ways
WHERE the_geom && ST_Buffer(ST_setSRID($1,4326), 0.1)
ORDER BY dist LIMIT 1 $$
LANGUAGE SQL;

Called like this:

SELECT * FROM next_way_point(ST_MakePoint(7.26, 47.52));

VS

SELECT gid, source, target, the_geom, ST_distance(the_geom, ST_SetSRID(ST_MakePoint(7.26, 47.52),4326)) as dist
FROM ways
WHERE the_geom && ST_Buffer(ST_setSRID(ST_MakePoint(7.26, 47.52),4326),0.1)
ORDER BY dist LIMIT 1;


texts:
I am trying to do point to point distances in QGIS but I have points to point distances that are being calculated across oceans.

Given people don't generally fly around cities does anyone know a way I can calculate the distance between points but restricted to land routes only (preferably roads would be great but that'd be an added luxury)?


texts: Simulating sea level rise using ArcGIS Desktop?

texts:
I am trying to look at how sea level rise would cover a certain area using a DEM.

Is there a way to do this sort of simulation using ArcGIS desktop?

I am trying to look at how much area would be covered with 1ft of sea level rise, 2ft, etc, so set heights.


texts: How can I automate interactions with vector features in OpenLayers3 (canvas) similar to Leaflet (SVG)?

texts:
I want to automate and test some aspects of my map application in OpenLayers 3. I mainly want to automate clicking on markers and other vector features(ideally through the DOM using XPath like SVG).
In OpenLayers 3 right now, this is almost impossible to do, since the vector features are not present in DOM due to the canvas. I can visually see them of course.

I tried doing this in Leaflet while using the SVG renderer for vectors and selenium web driver, and it was a piece of cake.

Is there any way to automate the canvas based vector layers right now?

Also would it be possible to render vector features in SVG instead of canvas in the future for OL3? (like it used to be in OL2)


texts: Script to join .csv file to shapefile and export changing column names

texts:
I have a shapefile with an attribute table like this:

FID   Shape      ID
0     Polgyon    3
1     Polgyon    3
2     Polygon    3
3     Polygon    3

and numerous csv files which look something like this:

FID   p-value  t-stat
0     0.50     0.0004
1     0.01     0.0006
2     0.02     0.0007
3     0.04     0.0008

and I am joining the csv files to the shapefile and then exporting like this:

import arcpy, os
arcpy.CheckOutExtension('Spatial')

arcpy.env.workspace=r'H:\path_to_csv_files'
shapefile=r'H:\pathway_to_shapefile.shp'
out=r'H:\output_pathway'
#loop through csv files
for table in arcpy.ListTables("*"):
#Create a new feature layer, this will remove previous join
  arcpy.MakeFeatureLayer_management(shapefile, "parcelLyr")
#   Join table to feature layer
  arcpy.AddJoin_management("parcelLyr", "FID", table, "FID", 'KEEP_COMMON')
  tableName = table.split(".")[0]
  if not os.path.isdir(out):
            os.mkdir(out)
#   #Export joined layer to new feature clas
  arcpy.FeatureClassToFeatureClass_conversion("parcelLyr", out, tableName)

but my out put table looks like this:

FID   Shape      ID  NDVI__5_cs   NDVI__5__1
0     Polgyon    3   0.50         0.0004
1     Polgyon    3   0.01         0.0006
2     Polygon    3   0.02         0.0007
3     Polygon    3   0.04         0.0008

NDVI_5.csv is the name of my first csv file so it is for some reason putting the name of the csv file as the column headers after the join.  When I do the join manually in arcgis this doesn't happen though.


texts: How to join CSV data (based on admin boundary names) to shapefile of aforementioned boundaries in R?

texts:
DATA:
CSV data,
Shapefile

Hi all!

I'm looking to perform what I've been taught is called a spatial join. After researching this further, it's looking like this operation is most often called a table join. So sorry for any confusion in the wording!

What I want to accomplish is very straightforward. I have a CSV that contains data for admin boundaries in the Philippines. I have a shapefile of said admin boundaries. I want to join these together in R to make choropleth maps. Easy enough, right?

To do this in QGIS, I simply add the shapefile and CSV as layers. Within the shapefile, Properties -> Joins, then add the join from there. Here's a screenshot for the sake of complete clarity:

As you can see, the Admin category of the CSV file matches the REGION category of the shapefile's attributes. Once those are joined, I am easily able to create choropleth maps from any of the other attributes from the CSV. Quick dummy-proof example here:

Now how the heck do I do this in R? It seems like it should be really easy to do, but I can't find much information, as all of my queries have yielded results of spatial joins concerning "join attributes by location" from data management tools in QGIS. I am not trying to find out which points fall within x polygons. I simply need to join the CSV to the shapefile so that I can make a choropleth map!

This question seems most similar to what I'm asking, but I couldn't sift through it enough to make sense of what I needed out of it.


texts: How to combine If Then with InStr in ArcMap field calculator?

texts:
I want to write an if then in field calculator that tests if a field in the attribute table contains one of several versions of another string (the versions vary by capitalization).

Trouble is, I can't get even a simple if then to work, let alone an if then else if, which i would need here.  I'm getting a 9999 error when I try the following:

(in the pre-logic script code)

Dim onlorain as Long
If InStr( [Export_Output.std],"Lorain" )>0 Then
onlorain = 1
else onlorain = 0
end if

(in the variable = ('parcel.onlorain =') box)

onlorain

I've posted a screenshot here:

I'd thought that ultimately, I'd be striving (when I build in the capitalization variants) for something like

Dim onlorain as Long
If InStr( [Export_Output.std],"Lorain" )>0 Then
onlorain = 1
elseif InStr( [Export_Output.std],"LORAIN" )>0 tHEN
onlorain = 1
elseif InStr( [Export_Output.std],"lorain" )>0 tHEN
onlorain = 1
else
onlorain = 0
end if

But if the 'simple', 1st version of this won't work, this won't work either.

Hope someone has thoughts.

Thanks in advance!


texts: How to eliminate points by difference (vector geo processing tool)?

texts:
I have a point that falls outside the bounds of my area of study. What I did was select it and resorted to vector commands - GIS - difference. Unfortunately I can not get the desired result. My shapefile has 500 points and I have to eliminate one and I can not get 499 . I have some doubts what to put in the vector layer of input and difference layer.

How do a remove this one point?


texts: What is Time of Validity in a GPS System?

texts: What is time of validity in a GPS system?  For defining time of validity, there are only a limited definitions available.  My interpretation for time of validity is the built-in system time (like on your computer) while GPS time is time from the GPS receiver.


texts: Buffer Polygon, BUT Inside and Only for a Distance

texts:
There has to be an easy way to do this. Been plugging away for an hour and I don't know if I am brain dead or what.

What I need to do is buffer a polygon on the inside, but only for 500 metres. If you have a look at this image:

You can see exactly what I don't want to do. The selected polygons (cyan) I buffered by -100 metres. What I want to do is actually create a buffer for the area from the polygon boundary to where it buffered it here.

I know I could do an erase on the original polygon with the above buffered result, but that seems like  an extra step and I want a one step solution.
I also know I could convert the polygon to a line and then buffer it, but once again extra step extra result.
Python programming could sure do this, but I really don't want to waste the time going down that road.

I do believe there is a way to do this, with one tool and one result, if not, well Esri should add that in.

I'm looking to know if there is an Esri ArcGIS Desktop tool to accomplish this?


texts: How to detect other objects are on which side of an object based on lat long

texts:
I am working on a system that detects the lat/long of vehicles. I want to track other vehicle's position: left, top-left, right, top-right, bottom-right, bottom-left etc based on lat/long.

I have used following logic to detect the bearing and direction of a vehicle:

https://www.dougv.com/2009/07/13/calculating-the-bearing-and-compass-rose-direction-between-two-latitude-longitude-coordinates-in-php/

http://www.movable-type.co.uk/scripts/latlong.html

Using the above script, I am able to detect the direction and bearing of other vehicles but could not get their position: left, right etc.


texts:
I have a list with over 10000 street names and places. These street names indicate where Radar Speed controls are located in Germany. I Need to have a map, where those Points are tagged. I did this for about 200 streets manually but it takes a lot of time.

Does anyone know a solution how to get this done automatically?

Does this work with QGIS and OpenStreetMap?


texts: Will ArcGIS 9.x tools work in ArcGIS Pro?

texts: I'm looking to use corridor designer, a tool for building wildlife corridors, but i has not been updates since the 9. something, while I have ArcGIS Pro 2, will the tool still work?


texts: Altering PostGIS table with X and Y values of lines

texts: I have a PostGIS table of a road network. The roads have topology and have been split with pgr_nodeNetwork. I also have source and target columns for each of the road segment. I want to be able to perform pgRouting with MapServer and OpenLayers. So I assume that the click in the OpenLayers WMS service will send coordinates to the server. Is there any way to alter the PostGIS table and ad x and y columns (like x1, y1 and x2,y2) without having to do it manually ?


texts:
I am trying to classify an image with SVM and RF in ArcMap. I have trained the classifiers and have the classifier definition files. I tried to classify but could not get RF and SVM from the drop-down menu for classification.

Can someone guide me how to classify with SVM and RF in ArcMap?


texts: WFS via Geoserver - Proxy or MapDocument Issues?

texts:
I struggle getting a WFS overlay from geoserver to work. I think I meet all the requirements for the WFS to be displayed in my map, but the request seems to get stuck somewhere... FireBug gives the error-message:

Access to restricted URI denied"  code: "1012
http://www.openlayers.org/api/OpenLayers.js
Line 1188

whereas Line 1188 contains the following:

this._object.send(vData);if(bGecko&&!this._async){this.readyState=cXMLHttpRequest.OPENED;
fSynchronizeValues(this)while(this.readyState<cXMLHttpRequest.DONE)
{this.readyState++;fReadyStateChange(this);if(this._aborted)

The whole applications are running on localhost, I installed Apache Webserver and configured a proxy since that was recommended in most Threads concerning WFS running on localhost. I used the proxy.cgi from http://trac.osgeo.org/openlayers/browser/trunk/openlayers/examples/proxy.cgi and allowed access for 'localhost','localhost:80' and 'localhost:8080'same with my IP (IP, IP:80, IP:8080).

When I visit localhost/cgi-bin/proxy.cgi I get to see the OpenLayers Page so that should be working!
My WFS Variable looks like this:

var wfs = new OpenLayers.Layer.Vector("WFS", {
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.WFS({
            version: "1.1.0",
            url: "http://localhost:8080/geoserver/wfs",
            featurePrefix: "workspacename",
            featureType: "layername",
            featureNS: "http://lvr", [tried many, it doesn´t need to be an existing URL right?]
            srsName: "EPSG:4326",
            geometryName: "the_geom"
            })
            });

then I set

OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url=";

Is that one complete and where do I put it in the map document?
Can anybody help? Tried a lot of options, nothing really helped but it seems like it´s some connection/security issue rather than smth related to the map script itself... Do I neccesarily need a format variable, a schema location or a stylemap? Anything else important missing?

Update from comments:

I can open it, typing localhost:8080/map.html... And when I just did that, the WFS data showed up!

Can´t believe it - I just accessed the html wrongly?! Beginners mistake. Thanks geographika...


texts: How can I identify regions with some points as centers?

texts:
I have a Postgres database of entire U.S. zipcodes and their lat/long data.

My requirement is to divide these zipcodes into non overlapping regions.

All regions have some center zipcode. All regions grow outward, starting from it's center until they intersect with another region.

EDIT: I do not have information about the radii/sizes of the regions. All I have is a list of zipcodes that can be considered as centers of the regions.

Is this achievable using Postgis?


texts: How to make transparent the Index 0 of the color table of a Tif/GTiff/GeoTiff using gdal?

texts:
With a Tif file, I'm trying to generate tiles using gdal tools.

gdalinfo

gdalinfo I_III_IVA_MRA_WGS.tif | more

Output

Driver: GTiff/GeoTIFF
Files: I_III_IVA_MRA_WGS.tif
Size is 17119, 23309
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (119.423193410175510,18.468686780527399)
Pixel Size = (0.000179663056824,-0.000188223140884)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 119.4231934,  18.4686868) (119d25'23.50"E, 18d28' 7.27"N)
Lower Left  ( 119.4231934,  14.0813936) (119d25'23.50"E, 14d 4'53.02"N)
Upper Right ( 122.4988453,  18.4686868) (122d29'55.84"E, 18d28' 7.27"N)
Lower Right ( 122.4988453,  14.0813936) (122d29'55.84"E, 14d 4'53.02"N)
Center      ( 120.9610193,  16.2750402) (120d57'39.67"E, 16d16'30.14"N)
Band 1 Block=17119x1 Type=Byte, ColorInterp=Palette
  NoData Value=0
  Color Table (RGB with 256 entries)
    0: 255,255,255,255
    1: 0,150,0,255
    2: 0,0,0,255
    3: 0,0,0,255
    4: 0,0,0,255
   ...
   22: 0,0,0,255
-- More  --

When I finally generate the tiles using gdal2tiles, it includes the canvass of the image.

My question is how can we make transparent the Index 0 of the color table of a Tif/GTiff/GeoTiff using gdal?


texts: Syntax error at or near "1" -- PostGIS in Action -- Chapter 11 SQL data dump


texts: If all my data is in UTM NAD83 Zone 10 and I get new data in GCS WGS 1984 do I have to reproject the WGS data to UTM NAD83 Zone 10 for it to align correctly or can I just bring it into an mxd that is already in UTM NAD83 Zone 10. If I need to reproject it using ArcCatalog how do I do this?


texts: I'm using qgis 2.16.  I am attempting to add additional data points (points only, no lines etc.).  The point symbol ( a simple circle) will not display, however the label (the value of the point) will display.  The added point data is added to the Attribute table.  How do I fix this so the points will display?


texts: Selecting distance-based features using ArcGIS ModelBuilder?

texts:
I have a shapefile file of rows and I need to make a selection based on distance, how to execute this select command automatically?

I need to select and export to a new layer the feature with greater distance.

I want to use ModelBuilder for this!

I use Arcgis 10.3.


texts: Removing overlapped buffers using Google Earth Engine

texts:
I have a line split up into line segments, with each segment having a buffer around that. I’m trying to remove where they overlap (basically where the darker circles are within the line buffer itself), so that when I later mask by buffer it doesn’t have duplicated results.
Anyone know how to remove those overlaps (and ideally keep the land they covered still within ONE buffered area)?
var lines = ee.FeatureCollection(
linetosplit.cutLines(ee.List.sequence(0, linetosplit.length(), 9000))
.geometries()
.map(function (geometry) { return ee.Feature(ee.Geometry(geometry));
}));

var buffered = lines.map(function (feature) {
return feature.buffer(2000);
 });

Map.addLayer(buffered);

print(buffered)

Link to code: https://code.earthengine.google.com/b6b989b9d0bed22a8882f95bd68d0d4d


texts: Define a projection for LAScatalog chunks using the lidR package in R

texts: Using lidR::grid_metrics with a LAScatalog is it possible to use a second projection for the processing "chunks" such that the resulting raster cells align with a raster image (sentinel 2)? I know there is the realignment option for a catalog to change the origin of the chunks, but it seems I need to take it a step further. I want to avoid re-projecting the imagery which would result in information loss. I could always re-project the LAS data, but this is undesireable as well. The EPSG of the LAS tiles is 2992 (Oregon Lambert Feet) and the imagery is 32610 (UTM 10N).


texts: Count number of group layers or Layers inside group layer using ArcPy?

texts:
I am looking for a way of counting the amount of layers within a group layer

I have an empty group layer inside my map and I am looking for a way of initiating a process if there are 0 layers inside the group layer.
So for I have been listing layers inside layer. But if there is nothing inside .. there is nothing to loop through.

Is there a simple way to count the amount of sub layers inside a group ?
For example I have a MXD with one empty group layer inside it. I want to count the amount of layers inside the group layer
So far I am only able to access the layer like this

import arcpy
group = []
mxd = arcpy.mapping.MapDocument(r"C:\Users\jackw_000\Documents\mxd\layer_mapV2.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
for lyr in arcpy.mapping.ListLayers(mxd, "my_group_layer", df):
     group.append(lyr)

print group

Returns ['My_group_layer']

When I print out the list it returns my_group_layer... . I am missing something.. ?
A sub note , I am running this from IDLE , so at no time I am in seeing whats inside the TOC interface


texts: Identify features of a group layer by GeoServer in QGIS

texts:
I encountered a strange behavior in QGIS version 3.8.1. A group layer consisting of raster, vector and PostGIS layers compiled via the GeoServer 2.15.1 is integrated as WMS into QGIS.
If I want to identify the features at a point using QGIS, only the information for the lowest layer is displayed in the Identify Results. In the layer preview of the GeoServer all features at this point are displayed.

Does anyone know if this is an error in QGIS or requires a specific setting in QGIS or the Geoserver?

Single layers are all queryable.


texts: Fast and robust open source portable script or tool for interpolation - triangulation or similar

texts:
I am looking for some fast and robust open source interpolation library/script/tool.

The goal is to re-interpolate relatively big (1M+ points) unstructured (xyz) grids into high resolution cartesian grids (tif files or other format I can process with GDAL) to produce nice smooth contours.

So far the only usable way I found is to use triangulation in saga_cmd but the more points I have the longer it takes and from some higher number it is not finished even after hours while if I do the same in other GIS software through GUI it takes only seconds to maximum few minutes.

I wanted to use GDAL but though nearest is fast, it doesn't get me detail I need and other methods took ages even for lower numbers of points... I am not sure if/how could I use GRASS. One of my requirements is that the solution must be portable (so far I used portable python, gdal and saga).

Based on answers so far:

Example of data

gdal_fillnodata is incredibly fast but leaving spikes no matter how much smoothed, same with SAGA close gaps, where the speed is highly dependent on tension threshold, the bigger the faster

all (2 and nearest) gdal interpolations are not really suited for this task - nearest just takes the nearest value (not really interpolation), moving average is one step better but creates big steps not suitable for contours, the IDW is very hard to set to some reasonable results and anyway it takes too long

so far the best results I got are still with triangulation. Bus as I mentioned it takes long time with SAGA for middle size grids and for big it fails

There are still few things I haven't tried yet or not tested properly:

SAGA Close Gaps with spline
python NUMPY and Skipy as mentioned in this thread:
Smoothing/interpolating raster in Python using GDAL?

Both I am not sure about the speed and will be glad for any experience or better sugestions. Thanks.


texts:
I am trying to follow the step-by-step directions in the answer to "how to separate Python installation with ArcGIS" but cannot find the path file, DTBGGP64.pth, in C:\Python27\ArcGIS10.2\Lib\site-packages.

Has someone has faced a similar issue?

I am wondering if there is a workaround by including a path file that tells Anaconda where arcpy is installed? I tried to create a .pth file containing the following text:

# .pth file for arcpy

C:\Program Files (x86)\ArcGIS\Desktop10.2\bin
C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy
C:\Program Files (x86)\ArcGIS\Desktop10.2\Scripts

and saved it to C:\Anaconda\Lib\site-packages (that is, if this is indeed where all the packages are saved) but I get the following error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\__init__.py", line 21, in <module>
    from arcpy.geoprocessing import gp
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\geoprocessing\__init__.py", line 14, in <module>
    from _base import *
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\geoprocessing\_base.py", line 14, in <module>
    import arcgisscripting
ImportError: DLL load failed: %1 is not a valid Win32 application.

I am running conda 3.7-py27 on a windows 7 (64bit) platform. Anaconda is installed in C:\Anaconda and ArcMap is installed in C:\ArcGIS\Desktop10.2. Is there an alternative solution to get arcpy to work with Anaconda? I just want to be able to run my scripts from Powershell.


texts: Accessing geonetwork API by JSESSIONID cookie

texts:
I am trying to execute this simple request on geonetwork, which should give me catalog statistics details (about private datasets for the corresponding user).

http://localhost:8080/geonetwork/srv/eng/q?_content_type=json&summaryOnly=true

It works, when I add basic authentication credentials, but it does not, when I add a JSESSIONID cookie.

I create a new Geonetwork session, getting this cookie back:

Set-Cookie: JSESSIONID=3038E369A5A499AF410F94ABF8A0B12D; Path=/geonetwork/

When I send this cookie with the above request, tomcat receives this request (from the catalina log):

"GET /geonetwork/srv/eng/q?_content_type=json&summaryOnly=true HTTP/1.1" 200 343 "serverTime=1489326597903; sessionExpiry=1489326597903; JSESSIONID=3038E369A5A499AF410F94ABF8A0B12D"

For me, this looks good. But geonetwork does answer the request just with the dataset statistics, that are public. Private datasets will be excluded from the response.

I examined the requests made in the Geonetwork UI and they look exactly like  what I am sending. Is there anything that I miss out at this point?

EDIT:

I am using version 3.0.3 of GeoNetwork.

The good point is: If I test from the Firefox RESTClient with a JSESSIONID created by logging in via the user interface, the session is actually recognized for private access.

So the problem seems to be in my code, where I create a new GeoNetwork session.

If I call with wrong username/ password, GeoNetwork gives me a 200 OK but no JSESSIONID.
If I call with correct username/ password, GeoNetwork will give me a 200 OK and a JSESSIONID. But apparently, the JSESSIONID does not seem to be qualified to access private datasets.

What I am basically doing to get the session id, is

Map<String, String> headers = new HashMap<>();
headers.put("Content-Type", "application/x-www-form-urlencoded");

String url = "http://localhost:8080/geonetwork/j_spring_security_check";
String body = "username="+userLogin+"&password="+userPass;

URLResponse r = HttpURLClient.sendRequest(url, headers, body, HttpMethod.POST);

Is that call to j_spring_security_check probably wrong?


texts: Problem Access-Control-Allow-Origin.GeoServer ISS OpenLayers

texts:
I use OpenLyaers in order to show france Departements. My website is generated in ASP .NET MVC 4 (i use Visual studio 2012) ans i use GeoServer to store WFS files.

I want to access to my WFS layer by

var wfs =  new OpenLayers.Layer.Vector("States", {
            strategies: [new OpenLayers.Strategy.BBOX()],
            protocol: new OpenLayers.Protocol.WFS({
                url: "http://localhost:8080/geoserver/wms",
                featureType: "states",
                featureNS: "http://www.openplans.org/topp"
            })
            });

map.addLayer(wfs);

In Google Chrome console i have :

XMLHttpRequest cannot load http://localhost:8080/geoserver/wms. Origin http://localhost:50431 is not allowed by Access-Control-Allow-Origin.

How i can parametre Geoserver to solve this problem ?


texts: Launching GRASS GIS 6.4.3 GUI on windows 8.1?

texts:
I downloaded GRASS GIS from this website：
GRASS GIS - Home, and I choose winGRASS 6.4.3 standalone installer to download.

After I installed, I click on its icon：

I can't launch it even though I have used 2 ways：

just left-click the mouse
right-click the icon and choose start as an administrator

can't launch means after I click the icon, the cmd icon shows up in less than one second on the windows bar below(even no windows come up), and then nothing comes up.

This really bothers me because I can't open GRASS!

I'm working under windows 8.1, is there someone using windows 8.1 facing the same problem like mine?


texts: Enabling DriveTime (Minutes) in Network Analyst Service Area Creation?

texts:
I'm using ArcGIS 10.1 Network Analyst Extension to try and do a drive time analysis for fire station service areas.

Using this tutorial I am going through step by step to make this happen, however, I am unable to create the analysis using drive time as an impedance and can only choose distance as an option:

My road network dataset has a field that shows the time it would take to traverse the segment at 40mph.

How do I enable DriveTime as an option in my analysis settings?


texts: Extracting values between neighbouring raster cells using QGIS?

texts:
I am trying to find out which green belts in a city are situated below (or almost on the same height as) the roads. I have a shapefile of all the paved roads and also a shapefile of the green areas. Also, I have a raster of 0,5m with the height values of the area.

First of all, it is important to know where the road and green belts adjacent. Therefore I did the following:

Created a buffer around the green
Created a buffer around the roads
Used intersect to find out where they overlap/adjecent

Next, I calculated the slope on the height raster and clipped this to the intersect vector layer that I just created. I now know how steep the slopes are on the areas where the green and roads meet. This is useful, but I think using a slope is not the best solution. This because I can now filter on areas where there is no to little slope which allows me to find the areas where green and road are almost on the same height. However, it does not really reveal the locations where the green belts are necessarily below the roads and also doesn't show by how many centimeter they are lower.

That is why I am trying to find a way to extract the height values of the cells where the green and the road adjacent. This might be a long shot, but if you know how I could solve this please let me know!

I use QGIS.


texts: Select duplicate attributes in same shapefile in QGIS

texts: I have this vector layer of lines between 2 nodes. I need to select all second entries that have the same length of the previous one. Can someone help me with this?


texts:
I am trying to toggle visibility of a GeoServer WMS layer.

The code below doesn't render the layer, regardless of if setVisibility is true or false

var countryParks = new ol.layer.Tile({
                 title: 'Country Parks',
                 source: new ol.source.TileWMS({
                 url: 'http://localhost:8080/geoserver/opengeo/wms',
                 params: {"LAYERS": 'opengeo:country_parks'},
                 serverType: 'geoserver'    }),
});
countryParks.setVisibility(true);
map.addLayer(countryParks);

And this code does render the layer, true or false regardless

var countryParks = new ol.layer.Tile({
                 title: 'Country Parks',
                 visiblity: false,
                 source: new ol.source.TileWMS({
                 url: 'http://localhost:8080/geoserver/opengeo/wms',
                 params: {"LAYERS": 'opengeo:country_parks'},
                 serverType: 'geoserver'    }),
});
map.addLayer(countryParks);

How does one toggle layer visibility, if not as in these attempts?


texts: How to reclass zero values for an aspect raster layer

texts:
reclass in QGIS 3.4.11 to reclassify a raster file of aspect values. All cells reclass apart from zero values which now show as 'no data'. Text file is as follows:

0 thru 23 = 2
23 thru 68 = 2
68 thru 113 = 1
113 thru 158 = 1
158 thru 203 = 4
203 thru 248 = 4
248 thru 293 = 3
293 thru 338 = 3
338 thru 360 = 2


texts: County Shapefile goes into body of Water

texts: I am trying to map three counties in QGIS with a layer of census tracts on top. The county border is in white and the census tracts are in green. Is there a way to not show the county border but leave the county names? It appears the census tract body of land matches what I see on California maps. It appears these three counties have borders that go into the water.


texts: QGIS "ERROR" written in all fields in attribute table of a View coming from PostgreSQL/PostGIS

texts:
I have a table of around 800000 records with POINT geometry in PostgreSQL/PostGIS. When I directly open the table in QGIS, the attribute table opens fine and the fields are populated.

But when I open the View of same exact table, "ERROR" is populated in all fields in the attribute table.

If I reduce the size of the table or reduce columns then the attribute table opens fine for a View as well.

I dont know what is the problem. Is there some memory limitations for a view in QGIS?


texts: Adding NetCDF-File to Mosaic Dataset in ArcGIS Desktop?

texts:
In ArcGIS Desktop 10.4 I have a NetCDF-File from this source http://eca.knmi.nl/download/ensembles/downloadchunks.php (grid 0.25, years 1995-2016), and am trying to add it to a Mosaic Dataset I already created. The data consists of temperature over time (in days).

I tried two variations:

reading the .nc-file as raster, saving as layer-file and then adding the layer-file to the Mosaic Dataset by the "add rasters to Mosaic Dataset"-tool.
reading the .nc file directly into the Mosaic Dataset using the "add rasters to Mosaic Dataset"-tool. Here I also ticked the variable I need (tg).

Both Variations gave the same result: the features are drawn, but the numbers are wrong. Don't know, where they are coming from, but it's obviously not taking the variable I ticked.

Does anyone have an idea what I did wrong?


texts:
I don't have any formal education in GIS so I'm thinking that I'm looking for something that exists, I just don't know what it's called.

I basically have a large raster data set that will get re-sampled into a pyramid of various resolutions (1m, 2m, 4m, 8m, etc) Is there a way to generate the pyramid with gdal using a custom projection?

And once that's created, is there a way to run operations on specific tiles and save them in a new pyramid in an organized way that can be referenced later?

I'll likely have multiple derivatives and will need to be able to check for the existence of a tile in a derivative pyramid, if it's not there, then create it based on the original pyramid. the tiles that get selected from the derivative pyramid would be based on a polygon overlapping an index for the pyramid of some sort.

I hope this makes sense. If this is something that I need to build from the ground up, then maybe someone can at least point me in the right direction so I know more about whats already available along these lines.

I should also mention that the derivative tiles will be made up of data of a dataset with tiles of different boundaries. (Such as combining a hillshade with aerial imagery) all the pyramids will all need to be built with tiles of exact same boundaries as the original dataset.


texts: gdal_translate with custom corner coordinates in python

texts:
I wish to convert this Python code using a gdal command line tool to using the gdal python library. I am using it to convert some PNG satellite images I am downloading to georeferenced TIFFs. I manually specify the bounding box of the images when I download them (in EPSG:3857), so I know all the info about the images, I just wish to store all the info in a single georeferenced TIFF file.

My current solution is to use gdal_translate called by subprocess but apparently this can all be done natively in python, so I would prefer to do that rather than rely on this hacked-together workaround.

from osgeo import gdal
import subprocess
import shlex

# ll_x etc... are all corner coords I calculate before hand
subprocess.run(shlex.split("gdal_translate -of GTiff -a_ullr {} {} {} {} -a_srs {} {} {}".format(
    ll_x, ur_y, ur_x, ll_y,
    "EPSG:3857",
    "input_sat_image.png",
    "referenced_sat.tif")))

I have seen individual questions that cover pieces of this, but I am struggling to put them all together.


texts: PyQGIS Raster Band Stats - Access pixel values of raster layer

texts:
I am trying to access individual pixel values in a raster layer I'm working with. I can use QgsRasterBandStats to get some information about the raster layer with the following code:

layer = iface.activeLayer()
print("Active Layer: ", layer.name())
provider = layer.dataProvider()
extent = layer.extent()
stats = provider.bandStatistics(1, QgsRasterBandStats.All)
print("min value =", stats.minimumValue)
print("max value =", stats.maximumValue)

I'm using the minimumValue and maximumValue attributes of QgsRasterBandStats. All the attributes are listed in the documentation.

In this raster layer I'm working with, there are really only 7 raster values. If I go into the symbology and manually classify these values, you can see them below.

What I'd really like to do is capture these values in a list, for example. Is there a way to access the individual pixel values in a raster object using QgsRasterBandStats?

Or could I choose a raster value (155 for example) and print the number of pixels with that value? Using something like the following pseudocode?

desired_value = 155
print(count of pixels with desired_value)


texts: Python error when running clipper on heatmap

texts:
I have created a heatmap from a delimited text layer and am trying to clip it to the boundaries of the map polygon. When I try to run the clipper I get this error:

  Traceback (most recent call last):
              File "C:/Users/Karongwe Base Manage/.qgis2/python/plugins\clipper\clipper.py", line 373, in run
                self.clip()
              File "C:/Users/Karongwe Base Manage/.qgis2/python/plugins\clipper\clipper.py", line 122, in clip
                features = layer.getFeatures()          AttributeError: 'QgsRasterLayer' object has no attribute 'getFeatures'
             2018-01-12T08:57:58    1   Traceback (most recent call last):
              File "C:/Users/Karongwe Base Manage/.qgis2/python/plugins\clipper\clipper.py", line 373, in run
                self.clip()
              File "C:/Users/Karongwe Base Manage/.qgis2/python/plugins\clipper\clipper.py", line 122, in clip
                features = layer.getFeatures()          AttributeError: 'QgsRasterLayer' object has no attribute 'getFeatures'

This is using a fresh install of QGIS 2.18.15. Can anyone help me? I am a novice with QGIS, I just know this process works on a different laptop here.


texts: Declarations required for Mapbasic api calls

texts:
I am trying to read a serial port by using api calls createfile, readfile, getcommstate, setcommstate etc...

this all requires the proper declaration of a structure "dcb"

typedef struct _DCB {
 DWORD DCBlength;
 DWORD BaudRate;
 DWORD fBinary  :1;
 DWORD fParity  :1;
 DWORD fOutxCtsFlow  :1;
 DWORD fOutxDsrFlow  :1;
 DWORD fDtrControl  :2;
 DWORD fDsrSensitivity  :1;
 DWORD fTXContinueOnXoff  :1;
 DWORD fOutX  :1;
 DWORD fInX  :1;
 DWORD fErrorChar  :1;
 DWORD fNull  :1;
 DWORD fRtsControl  :2;
 DWORD fAbortOnError  :1;
 DWORD fDummy2  :17;
 WORD  wReserved;
 WORD  XonLim;
 WORD  XoffLim;
 BYTE  ByteSize;
 BYTE  Parity;
 BYTE  StopBits;
 char  XonChar;
 char  XoffChar;
 char  ErrorChar;
 char  EofChar;
 char  EvtChar;
 WORD  wReserved1;
} DCB, *LPDCB;

My first question is how would you declare this in mapbasic?
Then how to determine the structures length?


texts:
In H2GIS, geometry can be converted to geography with function "ST_Transform"
For example:

SELECT ST_Transform(ST_GeomFromText(
'POINT(584173 2594514)', 27572), 4326);
-- Answer: POINT(2.1145411092971056 50.345602339855326)

How can I convert in the opposite way?
How can I convert lat/lng(2.1145411092971056 50.345602339855326) to a coordinate in Cartesian coordinate system?


texts: Digitized points are getting displaced automatically in QGIS

texts:
We are capturing 100's of drilling sample locations in QGIS. At times, we are copying the whole set of points from adjacent survey_lines to the next one. After that, we will capture the respective mineral values. But once we save it and open the shapefile again, the points are getting displaced randomly (but not all the points), just like in the image.

Actual values are 110,20,15,20,30 etc... But the respective points are displaced to some other place and had been replaced by some other points. One can see captured values (yellow) doesn't match with the actual values (black). Will it be a bug or due to any plugins or due to copying - pasting of points or any other thing?


texts: How can I loop through a list to create rasters and write them to .tif files?

texts:
I received 3 .mat files that each contain several daily maps. I need to take them and get them into ArcMap (10.2 if that matters) as a .tif with spatial information. I have already used the R.matlab package to read the .mat files and successfully did it for one day, however I do not want to have to write hundreds of lines of code to repeat it for the remaining days. Here is my successful code:

library(R.matlab)
library(raster)
library(sp)
library(rgdal)

setwd("C:/Users/Heidi Church/Documents/Thesis/FSLE")

FSLE2006mat <- readMat("FSLEheidi2006.mat")
FSLE2007mat <- readMat("FSLEheidi2007.mat")
FSLE2008mat <- readMat("FSLEheidi2008.mat")

#Test on single day in 2006
x <- raster(FSLE2006mat$F2006812, xmn=152, xmx=155.5, ymn=47.5, ymx=50.5, crs="+proj=longlat +datum=WGS84")
plot(x, xlab="Latitude", ylab="Longitude", xlim=c(152,155.5), ylim=c(47.5,50.5), main="Submesoscale Fronts: August 12, 2006")
writeRaster(x, filename="FSLE_20060812.tif", format="GTiff", overwrite=TRUE)

I had tried doing both a for loop and lapply, however kept getting this error:

Error in .local(x, ...) :
  unused arguments (xmn = 152, xmx = 155.5, ymn = 47.5, ymx = 50.5, crs = "+proj=longlat +datum=WGS84")

Here are examples of each that I tried:

for (i in FSLE2006mat[1:16])
{
FSLE2006[i] <- raster(x=FSLE2006mat[i], xmn=152, xmx=155.5, ymn=47.5, ymx=50.5, crs="+proj=longlat +datum=WGS84")
}

FSLE2006 <- lapply(FSLE2006mat, (raster(x, xmn=152, xmx=155.5, ymn=47.5, ymx=50.5, crs="+proj=longlat +datum=WGS84")))
FSLE2006 <- lapply(FSLE2006mat, raster)
FSLE2006 <- lapply(FSLE2006mat[1:16], raster)

I am wondering how I can loop through all of the elements in the .mat list for each year, make them a raster with the spatial info noted, and write it to a .tif file to import into ArcGIS. I should also note that 1:16 is to only use the first 16 of 18 elements, as the last 2 are not needed.


texts: How to enable anti-aliasing in ArcGIS for Desktop?

texts:
You can see from the screenshots that QGIS has, in my opinion, much more pleasant on-screen rendering of vector features.  I assume there is some sort of smoothing algorithm used in QGIS.  Is there a way to replicate this on-screen smoothing effect in ArcGIS?

(click images for full effect)

ArcGIS for Desktop ArcMap 10.2

QGIS 2.2 Valmiera


texts:
I want to hide WMS layer in Leaflet. But in
a condition, display again.
I tried removelayer. But once
It is removed, it can not add again.
here is the relevant code (generated with folium):
if (zoom > 9) {
  map_8e07456499fc4593a61b0bcabd1ab7cf.eachLayer(function (layer) {
    if (layer["options"]["attribution"] == "jeo1") {
      console.log("***   ", layer["options"]["attribution"]);

      var url =
        "http://localhost:8182/cgi-bin/mapserv?map=/map/jeo_yas2wms.map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=jeo_yas2&SRS=EPSG:4326&BBOX=" +
        bbox +
        "&WIDTH=1160&HEIGHT=540&FORMAT=image/png";

      map_8e07456499fc4593a61b0bcabd1ab7cf.removeLayer(layer);
    }
  });
}

I need normally get the map removed when I have zoom < 9.
What do you suggest?


texts: QGIS georeference raster to world file (WLD)

texts:
I am trying to georeference a raster (jpg). I can easily create a geotiff, but when trying to use the QGIS georeferencer to create a wld file instead, something goes wrong.
I have the target CRS set, and using liniar, the wld file is created.
But when I load into QGIS an error pops up and says the wld file contains no CRS, and thus the raster loads into the standard CRS (which is not the CRS I had to georeference the raster to). Anybody knows which knob I need to switch on or off? I can find plenty of pages explaining how to georef to geotiff, but this one I cannot find.

EDIT: I was following the guide here:
http://www.qgistutorials.com/en/docs/georeferencing_basics.html

first for the geotif, and then I tried to modify settings to save as a world file.
I set to liniar, set target crs to 4044 as pr the instruction on the page. But as was said, this is not saved to worldfile. and then I checked create world file. all filenames for .wld, .point and .jpg are identical.

using QGIS 2.14


texts: Finding average distance within polygon to point?

texts:
For my purposes, I can assume a homogeneous distribution of points within the polygon. What I want to find is the average distance from points within a polygon to a specified point. I don't know where to start and have found similar info online but nothing exact.

Ideally, I would like to do this for an entire Thiessen polygon layer. I have many polygons and a layer with one point in each polygon, and within each of those I would like to find the average distance to the point.

For example:

Within a country, what is the average distance one must travel to reach the Capitol? It can be assumed that people are evenly distributed throughout the country.


texts: Group features from feature layer by specific field without making multiple queries

texts:
I have a feature layer where one of the properties always takes 4 possible values.
Now I want to make 1 request to get all the features, group them by this property and feed each group to a different markercluster (or esri leaflet cluster) so I can have a separate cluster for every group. However it doesn't seem possible to do it without creating a cluster layer for each group explicitly, which would result in 4 separate requests and multiplying complexity. As far as I can see pointToLayer and iconCreateFunction doesn't really do the job (and they do weird things when points are removed or recreated), and panes only group the rendered markers, not the actual cluster itself.


texts:
I'm new to ArcObjects and I'm finding it difficult to work with the ArcObjects 'Feature classes'. Coming to my problem, I have a Map with few layers. I have coded for the layer names to display in the first Combo-box and the Column names of the Layer selected from the first Combo-box will be displayed in the second Combo-box. The values of the Columns will be displayed in the third Combo-box(So far, So good). But now I'm working on a query that displays the values lesser than the value selected in my third Combo-box. I have been using VB.Net coding for this purpose. The query code runs something like this.

Dim hook = ComboBox1.SelectedIndex(Table number)
    Dim FClass = ComboBox1.Text(Table Name)
    Dim Field = ComboBox2.Text(Column Name of the table selected)
    Dim Value = ComboBox4.Text(values in the Column selected)
    Dim pMxDoc As IMxDocument = TryCast(m_application.Document, IMxDocument)
    Dim pMap As IMap = pMxDoc.FocusMap
    Dim pLayer As IFeatureLayer
    pLayer = pMap.Layer(hook)
    Dim pFeatureClass As IFeatureClass = pLayer.FeatureClass
    Dim pTable As ITable = pFeatureClass

    Dim FilePath As String
    FilePath = "D:\Harsha\ArcMap_CustomTools\EGISDATA\"
    Dim pWSFact As IWorkspaceFactory
    pWSFact = New ShapefileWorkspaceFactory
    Dim pFeatureWorkspace As IFeatureWorkspace
    pFeatureWorkspace = pWSFact.OpenFromFile(FilePath, hook)
    pLayer.Name = FClass
    pTable = pLayer
    Dim pQFilter As IQueryFilter
    pQFilter = New QueryFilter
    pQFilter.SubFields = Field
    pQFilter.WhereClause = "'" + Field + "' < '" + Value + "'"
    Dim pSelSet As ISelectionSet
    pSelSet = pTable.Select(pQFilter, esriSelectionType.esriSelectionTypeHybrid, esriSelectionOption.esriSelectionOptionNormal, pFeatureWorkspace)
    pMxDoc.ActiveView.Refresh()
    Dim pFeatureCursor As IFeatureCursor
    pSelSet.Search(pQFilter, False, pFeatureCursor)
    Dim iRecCount As Integer
    s = vbNullString
    Dim pFeature As IFeature
    pFeature = pFeatureCursor.NextFeature
    While Not pFeature Is Nothing
        iRecCount = iRecCount + 1
        s = s & vbCrLf & pFeature.Value(Value)
        MsgBox(s)
        pFeature = pFeatureCursor.NextFeature
    End While

The problem is pFeature is always Nothing and it does not enter the while loop


texts:
This is a General Question on a possible workflow for Survey Data directly to ArcGIS

In most GIS/Survey workflows there is a field collection team and backoffice GIS team.

When the availability of high-end mapping devices are there, why can we directly edit a Geodatabase in a Survey device or rather create a GIS feature in Survey device itself, eliminating the CAD and GIS feature creation ?

What my Question is , is there any Survey Device ( Totalstation etc) currently supporting data  editing in GIS format ? Is there any Survey device capable to support a full fledged Geodatabase ?


texts: Does the inputs extent have to be identical for raster::intersect to crop a RGB GEOTiff?

texts:
Background: I have a shapefile with multiple polygons over a large area and then I also have GEOTiff image files. I want to crop the TIFF files with each of the polygons. The problem is that when I use raster::intersect to try and do it, I get the error:

Error in .local(x, y, ...) : extents do not overlap

This is my (pseudo)code for loading the shapefile and selecting a polygon from that shapefile:

 Mineralrights.shape<-readOGR(dsn = "C:/Users/blablabla/Shapefiles",layer = "Mineral Rights")
 Mineralrights.shape <- spTransform(Mineralrights.shape, "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")
 tempshape <- Mineralrights.shape[Mineralrights.shape$GID == "2097674",]

And this is my (pseudo)code for loading the GEOTiff as a raster:

  Aerialfiletocheck <- raster(paste0(getwd(),"/","2923AD_06_2014_783.tif"))
  projection(Aerialfiletocheck) <- projection(tempshape)

I don't think the projections is the issue here, since all the source files (shapefile and TIFFS) uses a standard WGS84 ellipse, I just set the projections explicitly for safety. Then call:

  clip <- raster::intersect(Aerialfiletocheck,tempshape)

And this is when I get the error about the extents that do not overlap. To be sure, I tried to simplify the GEOTiff to just its extent and then plotting that over the specific polygon I have selected:

  raam <- projectExtent(Aerialfiletocheck,"+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")
  raam.shape <- as(raam@extent, 'SpatialPolygons')
  plot(raam.shape)
  plot(tempshape,add=T)

Yielding this plot, where the large rectangle is the extent of the GEOTiff and the funny triangle being the polygon I selected.

These are the parameters requested

 > Aerialfiletocheck@extent
    class       : Extent
    xmin        : 24127.25
    xmax        : 29297.25
    ymin        : -3248250
    ymax        : -3242397
 > tempshape@bbox
            min       max
    x  23.25377  23.25875
    y -29.33836 -29.33380

Problem: To me it looks like their extents very much overlaps, hence my question: Does the inputs' extent have to be identical for raster::intersect to crop a RGB GEOTiff?
That would make no sense to me, because the whole idea of raster::intersect is to "crop" the one with the other, which would imply that their extents will overlap, but not necessarily be identical.

What else have I tried:
1) I have replaced the raster with another polygon and then get the expected result.

2) I have tried base::cropwhich give the same error, which I expected because I think when one input is a raster, intersect just calls crop

3) I have tried rgeos:gIntersection which results in an error about the raster not having a proj4string slot, which makes no sense to me.

4) I have read these posts (and many others): Extracting intersection areas in R

Crop a raster file in R

Extract Raster from Raster using Polygon shapefile in R

GDAL: Convert specific RGB (GeoTIFF) color to transparent


texts: Is it possible to get an entire street via a single Overpass API query?

texts:
I'm working on a game (http://igotchareverse.com/) where you build buildings along real world roads. Currently investigating a feature where a player who controls the most buildings along a street gets points based on the streets length.

The problem i'm having is that streets as humans understand them do not exist as data in OSM. The best way to get a street is to traverse the ways with the same names that are connected by nodes, taking some wild cards like roundabouts into account.

The main question is, can i get an entire street via OverpassAPI with a single query using one of it's ways osm_id or similar value? A street in this case is a series of ways with the same name connected by nodes.

Things that need to be considered are that many ways have the same name where we work, Sweden. Almost every major city has a Kungsgatan for example (Kings Street), so it's not simply a matter of getting a list of all the ways with the same name and calling it a street, they need to be connected with nodes.

Another issue is that there is to my knowledge any other reference or value that binds these ways together.

I can without a problem build function in my server that achieves this using multiple queries, but that seems not performance efficient and also i'd like to avoid making my own endpoints and changes to Overpass, if possible. Performance is also the reason i won't be making multiple calls like that from the client which will be running on a mobile device.


texts: I am working on a Python package that has a mixture of command line tools and ArcToolbox tools.  For Python development it is really useful to use VirtualEnv, especially when working with installers.  Is there a way to have either ArcMap or ArcCatalog recognize a current VirtualEnv session in Python, so that when a custom tools is executed it is using the VirtualEnv Python paths and packages instead of the default Python paths and packages at c:\Python26\ArcGIS10.0?


texts: How to style specific polygons from a geojson with Leaflet

texts:
I have a World Countries geojson, that has country codes as "ID" (feature.id) and names as feature.properties.name

I would like to color specific countries for making a "countries visited" map.

What would be the easiest way to provide a list of ids or names to assign a color to them?

Currently my code for one country is this:

 var map = L.map('map').setView([30.45,24.61], 3);
 function style(feature) {
        if ("Latvia" == feature.properties.name) {
            return {
                weight: 2,
                opacity: 1,
                color: 'white',
                dashArray: '3',
                fillOpacity: 0.7,
                fillColor: '#FFA000'
            };
        } else {
            return {
                weight: 2,
                opacity: 1,
                color: 'white',
                dashArray: '3',
                fillOpacity: 0.7,
                fillColor: '#000'
            };
        }
    }

var geojsonLayer = new L.GeoJSON.AJAX("countries.geo.json", {
            style: style
        });

geojsonLayer.addTo(map);

I would like to replace the "IF" thing with some sort of comma separated list of ids or names.

Here is a sample of that geojson file:

{"type":"FeatureCollection","features":[
{"type":"Feature","id":"AFG","properties":{"name":"Afghanistan"},"geometry":{"type":"Polygon","coordinates":


texts:
I added an image, that explains my demand (see Schiefelsbach river).
How can I do this?


texts:
I have installed qgis 2.0, I need to install add-ons and NEXT I get the error 'The supplement package decompression failed. probably is broken or missing in the repository. might also want to make sure you have write permission to the plugin directory:
/ home/cordero/.qgis2 / / python / plugins "I searched the file path indicated by the error but apparently I have not the directory, and uninstall and reinstall qgis on my pc but the problem persists


texts: Receiving Albedo values from Landsat with Google Earth Engine

texts:
I am new to Google Earth Engine. I would like to extract the Albedo for certain areas of Greenland. I am doing it with MODIS as well as Landsat and was wondering the following about Landsat:
I found that you use the following equation for albedo:
((0.356*B1) + (0.130*B2) + (0.373*B3) + (0.085*B4) + (0.072*B5) -0.018) / 1.016

Does it make more sense to use 'USGS Landsat 8 Collection 1 Tier 1 TOA Reflectance' or 'USGS Landsat 8 Surface Reflectance Tier 1' for that?


texts: How to update the layers in Geonode with cascading WMS from Geoserver?

texts:
We began a project in my office in which we would like to have all our geographical data in a central repository where we can search them, insert and edit metadata and view them on a map, so our choice is Geonode.

We have a lot of data around on different machines and servers and on different physical locations so we dont necessarily want to put them on a single drive and / or convert them into PostGIS layers. So what we were thinking was to have a couple geoserver instances locally on the different offices and then have geonode consume WMS services from these geoservers.

So what we did was to create cascading WNS on the geoserver instance that runs Geonode.
The problem is that we cannot get Geonode to read those new WMS services. We tried to run geonode updatelayers. The script runs without any error but only updates the existing shapefile layers and not the cascading WMS. When we run

 geonode updatelayers -s Ktimanet2

and the name of the store it gives an error:

 FailedRequestError: No store found named: Ktimanet2

This is what my layers look like in Geoserver:

Any ideas anyone??


texts:
I'm trying to draw polyline to connect an array of waypoints. But I'm getting an empty div like this

instead of a map with route.

My code is this:

where I've made a mistake?

Here I get an array of latitudes. The trouble is in these lines:

L.Routing.control({
            polyline: wps,
            lineOptions: {styles: [{color: 'black', opacity: 0.15, weight: 9},
                {color: 'white', opacity: 0.8, weight: 6},
                {color: 'blue', opacity: 0.3, weight: 15}],
                addWaypoints: false},
            createMarker: function() { return null; },
            show: false
        }).addTo(map);

If I write waypoints instead of polyline it works, but I need to make it with polyline.


texts: How to get OpenLayers map object in HTML form?

texts:
There is a form to within a web page, map of OpenLayers is displayed in the form. I look at the source code of Html, and further a hierarchical state in the form. My goal is to get the Polygon at the bottom, but do not know how to get the map object of OpenLayers.

Html source code

<form id="map_form" method="get" action="#">

          <svg id="OpenLayers.Layer.Vector_67_svgRoot">
           <g id="OpenLayers.Layer.Vector_67_root" style="visibility: visible;" transform="">
            <g id="OpenLayers.Layer.Vector_67_vroot">
             <path id="OpenLayers.Geometry.Polygon_3503" d=" M ...
             ...

Also tried the following methods.

1. map = new OpenLayers.Map (document.getElementById ('OpenLayers.Map_7_OpenLayers_ViewPort'))
2. map = new OpenLayers.Map (document.getElementById ('map_content'))
3. map = new OpenLayers.Map (document.getElementById ('map_object'))

Which method but also map.layers.length was 0.

How can I get the OpenLayers's map object?