texts: ArcPy RasterUniqueValueColorizer not working correctly

texts:
I have a raster layer which consists of unique values 0-4 and I want to color them through ArcPy.
I have tried this:
import arcpy
import os
arcpy.env.overwriteOutput = True
from arcpy.sa import *

#set the pathway to the maps
aprx = arcpy.mp.ArcGISProject(rel_path)

#get the map of interest
one = aprx.listMaps()[0]

#add new layer to map
one.addDataFromPath(new_layer_location)

#get all layers
one_layers = one.listLayers()

#get the first layer
l = one_layers[0]

#build attribute table so I can colorize it
arcpy.BuildRasterAttributeTable_management(l, "Overwrite")

#get symbology
sym = l.symbology

#make a dictionary for changing colors
legend_dict = {'0': {'RGB': [255, 85, 0, 0]}, '1': {'RGB': [0, 92, 230, 0]}, '2':{'RGB': [132, 0, 168, 0]}, '3':{'RGB': [56, 168, 0, 0]}, '4':{'RGB': [255, 255, 255, 0]}}

#set the colors to use based on value
if hasattr(sym, 'colorizer'):

    if sym.colorizer.type == 'RasterUniqueValueColorizer':

        sym.colorizer.field = 'Value'

        for grp in sym.colorizer.groups:

             for itm in grp.items:

                #change the color
                itm.color = legend_dict[itm.label]

#update symbology
l.symbology = sym

#save output
aprx.saveACopy(relpath + r'\\SavedOutput.aprx')

but when I look at this output all of the colors in the legend are transparent are have the information of:
{'HSV': [0, 0, 100, 0]}

I used this information to get this far:
https://pro.arcgis.com/en/pro-app/latest/arcpy/mapping/rasteruniquevaluecolorizer-class.htm
but I am not sure where it is different than what I am doing.


texts: QGIS - KML File Not Displaying Correctly After Conversion

texts:
I have a ESRI shapefile made up of 140 polygons that I want to display on Google My Maps. I have my project and layer CRS set to WGS84 which is standard for displaying spatial data on google. However, my polygons do not appear when I open the file on QGIS and when I try to upload it to My Maps they're displayed as random points/lines instead of polygons and are far off the map from their original location. My attribute table contains only the names of each individual polygon, some have the same name. What should I do to have my .kml display properly? I thought it would contain all the relevant metadata already.

-


texts:
A client of my company has asked us to add a CAD file to a WMS. We have a subcontracted IT company that will adress the WMS part, but we need to provide them with a geoJSON of our client's CAD file.

I have managed to transform the provided .DXF file into a shapefile but it seems I am not able to convert that shapefile into a raster properly. I need to do this step in order to georeference the raster and then export it as a geoJSON.

The problem is the following (please understand that I'm not a programmer and I have little to no knowledge in linux commands), using GDAL tools within QGIS, specifically gdal_rasterize. I give the program an input file, create a new output raster with resolution 1,00 horizontal and 1,00 vertical which gives me the following command:

gdal_rasterize -tr 1.0 1.0 -l 70080-00-ytg-mda-tr-001 <input_filename> <output_filename>

First of all QGIS throws an error telling me that I must input an extra parameter, shall it be -burn, -a or -3d. I don't know how any of this commands work, but after a little online research, I add -burn 0 so the code becomes the following:

gdal_rasterize -burn 0 -tr 1.0 1.0 -l 70080-00-ytg-mda-tr-001 <input_filename> <output_filename>

Then, the disk's free space is not enough (which is not true). Again, a brief online research tells me to add this at the end of the command:

--config CHECK_DISK_FREE_SPACE NO

Then the process finally starts doing something, but after a few seconds it stops and a new error appears: "The process was interrupted a while after starting correctly" (may not be the exact words, as I am a spanish user, I translated it).

I am given an output file but it is corrupted, and at this point I don't know how to continue or how to fix the error. I hope there is someone kind enough to give me a solution or guide me through the learning process of doing so myself.


texts:
I am trying to clip a raster with a vector layer mask (shapefile) using the clipper tool (gdal) but I've got the next error:
Traceback (most recent call last):

          File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\GdalTools\tools\widgetPluginBase.py", line 88, in onFinished
            self.base.onFinished(exitCode, status)
          File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\GdalTools\tools\dialogBase.py", line 208, in onFinished
            msg = unicode(self.process.readAllStandardError())
        UnicodeDecodeError: 'ascii' codec can't decode byte 0xf1 in position 85: ordinal not in range(128)

Both layers are in the same coordinate reference system, so this is not the problem


texts: Water reflectance effect on Sentinel 2 images

texts:
I'd like to know what is causing the effect on the water depicted on the following screen (sentinel tile id: S2B_MSIL1C_20190805T144739_N0208_R139_T20PPA_20190805T162922 located near Saint Lucia).
You can see a darker area on water were light seems to be decomposed resulting in heap of pixels of different colors (green, red, blue).
I've checked the data quality report but this effect was not mentioned.


texts:
I have a set of geotif images that I want shifted by an integer number of pixels in X and Y. However, those images are quite large, so I don't want to duplicate them (and I need to keep the original unmodified), so modifying the geotransform is not a solution for me.

Is it possible to create a set of vrts (with different shift values) without modifying or duplicating the original images?


texts: Calculating Mercator coordinates from lat/lon

texts:
I am new to GIS but I am working on a project mapping a building using floor plans and JOSM.

Using the PicLayer plugin I added the .png floor plan as a layer.

I wrote a .cal file to set the scale, rotation, and center position of the image.

In the JOSM default Mercator projection, the image moved only slightly when the center was changed (using lat and lon decimal degrees). The image was still centered somewhere around zero. In WGS84 the image moved to the correct location but was sheared and stretched.

I found this reference on converting lat/lon to a Mercator projection: http://wiki.openstreetmap.org/wiki/Mercator

The Python code under Elliptical Mercator correctly converts my longitude, but the latitude is off by about 0.5% (a few metres).

import math

def merc_x(lon):
  r_major=6378137.000
  return r_major*math.radians(lon)

def merc_y(lat):
  if lat>89.5:lat=89.5
  if lat<-89.5:lat=-89.5
  r_major=6378137.000
  r_minor=6356752.3142
  temp=r_minor/r_major
  eccent=math.sqrt(1-temp**2)
  phi=math.radians(lat)
  sinphi=math.sin(phi)
  con=eccent*sinphi
  com=eccent/2
  con=((1.0-con)/(1.0+con))**com
  ts=math.tan((math.pi/2-phi)/2)/con
  y=0-r_major*math.log(ts)
  return y

My nodes and ways were lined up correctly with the image with this latitude and longitude in WGS84, so I think they are correct and the conversion is not. I am in the northern hemisphere and using an area near 49 N, -122 E.


texts: Google Earth Engine: how to export a sloping image without zero patch

texts:
I want to export a strip image in GEE, and the image I want to export looks like:

The export code is:

    Export.image.toDrive({
      image: img,
      description: descp,
      scale: 30,
      folder: folder,
      fileFormat: 'TFRecord',
      region: img.geometry(),
      formatOptions: {
        'patchDimensions': [256, 256],
        maxFileSize: 104857600,
        compressed: true,
      },
      skipEmptyTiles: true,
    });

However, when I read the exported tfrecord files, I still get a lot of image patches which are full of zero values. So I'm wondering how can I avoid this kind of situation?


texts: Changing hasNoDataValue from TRUE to FALSE by using rgdal

texts:
I was wondering if there is a way in R to change hasNoDataValue from TRUE to FALSE.
Basically I have a raster that, when using rgdal::GDALinfo('myraster.tif'), I can see that hasNoDataValue" = FALSE:

But after some analyses with the raster package in R, when I export the output raster, the hasNoDataValue is automatically set to TRUE, and I need it to be FALSE, but I don't know how to change this.
Is there any way to do this in R?


texts: Stamp out and transfer values from one overlayed vector layer to another

texts:
I have a vector layer with administrative boundaries and another vector layer with postal code boundaries, which are more precise/smaller entities/objects.  I now want to transmit the values of one column/attribute of the first layer to the second.

Is it possible to overlay them and kind of stamp out this attribute and thus transfer it to the postal code layer?


texts: Filtering a polygon by area with multiple iterations?

texts:
I am currently working on a habitat conservation map in ArcGIS 10.3.  I want to filter out smaller areas from my polygon layer, but as I'm not sure what the appropriate habitat patch size is for this, I plan to perform multiple filtering operations.

For example, I'd like to filter out patches that are less than 5 acres in area and yield a layer output; then I'd like to repeat this filtering process with layer outputs for 10 acres, 20, 50, and 100.  I can do this manually, but it seems tedious enough that I could set this up in ModelBuilder or Python-but I need help with setting this up, as I am unsure of how multiple iterators will work in ModelBuilder with ArcGIS 10.3.


texts: High Resolution Monthly or Annual Precipitation data for four corner states or the Navajo Nation

texts:
I was wondering if anyone on here might know where I can find 30 meter resolution Average Annual or Average Monthly Precipitation data, that would cover the 4 corners area (UT, AZ, NM, CO) or basically the Navajo Nation?

I am working on a Mule Deer Habitat Suitability analysis for one of my class project's and have elevation and landcover data that is 30m resolution.

I wanted to try to do my analysis with 30m precipitation data too, from at least two different years or timeframes (i.e. 2008-2011 and 2014-2017 -OR- 2010 and 2013).

I know PRISM data is at 800m resolution, and could resize the cells to 30m, but would be great if I could use 30m precip data.

I've been looking around for a couple days now.


texts: ArcGIS - viewshed calculation using a DSM and DEM

texts:
I am fairly new at using ArcGIS.  For work we use it only to create viewshed analysis using DEM - or bare earth terrain.

I can create these fine, however I need to create one showing the effects of vegetation and buildings. I have received a DSM of the site, but if I take it through the normal process, as we do for the bare earth, it calculates the DSM surface, rather than calculating the view from the terrain.

So for example, where there are clearly trees and buildings, it shows that you can in fact see the observation point, even though you know you wouldn't because of the tree and buildings would be covering it.

I can't seem to find any step by step tutorials that can help me out,  I'm sorry for my lack of detail and wording on explaining my problem.
I am really really new to this, and have been thrown in the deep end by work wanting me to find the answer!


texts:
I'd like to create a map of sinkholes in ArcGis 10, I have a DEM, 2m resolution, is there a tool to detect local minimum points ? tried "sink" tool without any success, result is very small holes.

And I did a fill to a certain depth, and substracted it frome the original DEM, it worked but not as precise as I want and lots of artefacts, very small holes, so I need to find a way to define to what extent (diameter, depth) holes need to be filled and I did not find that for now. can you help me ?


texts: How to Display Coordinates of an Existing Point Layer in ArcMap

texts:
I have been asked to provide a list of coordinates for a point layer in our GIS.  I am currently using a state plane coordinate system.  How do I view/list this data.  Below is the data request in it's original form: (Each serial number is assigned to a specific city address)

If I were to provide you with a list of ONT serial numbers like this one:

Serial Number Lon Lat
S050731154
S031620362
S022960010
S030560404
S021860376
S043371128
S032140857
S030490038
S023080296
S022380211
S030150116
S023080386
S043371132

Would you be able to fill in the longitude and latitude of the ONT owning that serial number based on the data stored within the GIS?

Or on the other hand, do you have the ability to export such a list?


texts: Adding value to a polygon layer with a feature of the nearest polygon of another layer

texts:
I have two polygon vector layers, shown below in purple and orange. In the orange layer, each feature has a height value that I want added to the nearest purple feature.

Ideally, a new vector layer is made consisting of the features from the purple layer with the added height value from the nearest feature from the orange layer. But updating the purple layer will also suffice.


texts: Using Field Calculator vb to check for name with /n in the string?

texts:
I am calculating values for a field (Score) based upon the names from another field (Modified). One name is "blaugh\nblaugh". Field Calculator keeps treating the name as having a new line command.

I have tried using Like, Right and Contains but the Field calculator does not seem to be able to use those commands. I am stumped.

Ex1

elseif  [Modified] ="blaugh\nblaugh" then
Score =7

Ex2

elseif Right( [Modified], -6) = "blaugh" then
Score =7

Ex3

elseif [Modified] LIKE "*nblagh" then
Score =7

I am using vb in ModelBuilder and ArcCatalog 10.3.1


texts: Calculating the outer boundary of several Geometry objects in SQL Server 2008

texts:
I have lots of Polygons of datatype Geometry in SQL Server 2008. The image below shows how a select of all these Geometrys looks visualized.

What I need to do is create a Polygon which represents the outer boundary of all these polygons. So I used the response given to a previous spatial question I asked to create the following code:

DECLARE @test TABLE(geom GEOMETRY);

INSERT INTO @test SELECT geom FROM ForceBoundary

DECLARE @geom GEOMETRY
SELECT @geom = (SELECT TOP 1 geom FROM @test)
SELECT @geom = @geom.STUnion(geom) FROM @test

SELECT @geom

This produced the following result, which has cracks in it due to holes between the polygons:

So I updated my query with the following change:

INSERT INTO @test SELECT geom.Reduce(0.001).STBuffer(100) FROM ForceBoundary

Which improved the result, however it's not solving the issue completely and it also damages the outer boundary accuracy.

What is the correct way to achieve this? From looking through a list of the STxxxx functions I couldn't see one that seemed to provide the results I need?

I've tried SELECT @geom.STExteriorRing() but it doesn't return anything. Appending it with .STAsText() returns null.


texts: Map does not appear after following ESRI example of flex api

texts:
I followed the example from ArcGIS resource center website: https://developers.arcgis.com/flex/ to create a map, unfortunately after I run the flex builder application, the map does not appear. The following is my code and what happens after running the application. What's wrong with my project? I have tried to use IE, Chrome and Firefox to open the url, but the result is the same:

<?xml version="1.0" encoding="utf-8"?>
<s:Application
xmlns:fx="http://ns.adobe.com/mxml/2009"
xmlns:s="library://ns.adobe.com/flex/spark"
xmlns:mx="library://ns.adobe.com/flex/halo"
xmlns:esri="http://www.esri.com/2008/ags">

<esri:Map level="3" wrapAround180="true">
    <esri:center>
        <esri:WebMercatorMapPoint lon="0" lat="0"/>
    </esri:center>
    <esri:ArcGISTiledMapServiceLayer
        url="http://server.arcgisonline.com/ArcGIS/rest/services/
        World_Street_Map/MapServer"/>
    <esri:ArcGISDynamicMapServiceLayer
        url="http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/
        Demographics/ESRI_Population_World/MapServer"
        alpha=".70"/>
 </esri:Map>

</s:Application>


texts:
I like to work with the label toolbar in qgis, because it is intuitive and I can place labels manually.
However, I lose the position and rotation of labels very often in my cartography workflow. That means I have to repeat the whol process of positioning labels again.
What could be the reason behind that? And how can I prevent it in my future qgis map design workflows?


texts:
The Problem
I have a SVG of a triangle that I would like to fill my polygons with. The problem of just using the SVG Fill tool is the triangle would not overlap fill in the white space as shown bellow.

My Attempt
Therefore, I have attempted to use the  and  tags to repeat the pattern in the square SVG canvas, something like this:

<svg xmlns="http://www.w3.org/2000/svg" width="288.075" height="249.628" viewBox="0 0 288.075 249.628">
  <defs>
    <pattern id="kumiko" viewBox="0,0,288.075,249.628" width="100%" height="100%">
        <path d="M.037 249.12 144.037 0l.293.212L288.075 249.48l-.33.148L0 249.48l.037-.36 144-82.8L144.33.212l-.293 166.108 143.707 83.308-143.707-83.308-35.577-20.584 35.577 20.584 35.615-20.519-35.615 20.519-.037 41.103.037-41.103" fill="none" stroke="#000" stroke-width=".72"/>
    </pattern>
  </defs>
  <rect width="288.075" height="249.628" fill="url(#kumiko)" stroke="url(#kumiko)"/>
  <rect width="288.075" height="249.628" fill="url(#kumiko)" stroke="url(#kumiko)" transform="rotate(180, 144.0375, 124.814) translate(144.0375, 0)"/>
  <rect width="288.075" height="249.628" fill="url(#kumiko)" stroke="url(#kumiko)" transform="rotate(180, 144.0375, 124.814) translate(-144.0375, 0)"/>
</svg>

I previewed it with plugins in VS Code and it looks alright. However, when I apply this SVG into QGIS, it only resulted in polygons with no fills at all.

Does QGIS support SVGs with the pattern tag or is this a mistake in my code?


texts: Properly compressing KML file in Linux OS

texts:
I am trying to compress a KML file using gzip on a Linux system

gzip -cvf somefile.kml > somefile.kmz

When I go to try and verify the new KMZ file in Google Earth, I receive the following error:

  Open of file "somefile.kmz" failed: Parse error at line 1, column 0: not well-formed (invalid token)

What happened?


texts:
I am working on some scripts with 3D polygons using ArcPy (ArcGIS 10). Reading and data manipulation works fine.  However, writing the output back into a shapefile does not seem to work, it only generated 2D files. For writing I use the following function:

import arcpy
def write_new_MP(data,outfile):
        array = arcpy.Array()
        point = arcpy.Point()
        features=[]
        for feature in data:
            for part in feature:
                point.X=part[0]
                point.Y=part[1]
                point.Z=part[2]
                array.add(point)
            MP=arcpy.Polygon(array,None,True,True)
            array.removeAll()
            features.append(MP)
        arcpy.CopyFeatures_management(features, outfile,"", "0", "0", "0")
        return

the data is a list likedata=[[[5,5,5],[10,5,5],[10,10,5],[5,10,5],[5,5,5]],[[15,15,15],[10,15,15],[10,10,15],[15,10,15],[15,15,15]]]

Does somebody know why the Z-values aren't used when generating the output?


texts: What GIS analysis techniques were used in determining Gale Crater as the landing site for the Curiosity Rover?

texts:
What specific GIS analysis was used in determining Gale Crater as the landing site for the Curiosity Rover?

I'm assuming factors such as slope, elevation, ruggedness, etc were used as well as potential for finding water (what determines this).

I'm looking for a journal article (if it exists) or some technical document outlining the processes used for choosing Gale crater over other potential landing sites.

This is just for my own curiosity (pun intended).


texts: Combining vector and raster tiles in mbtiles file?

texts:
I have 4 mbtiles files..

OSM Streets for the State of Ohio (vector)
Hillshading for the State of Ohio (raster)
Contours for the State of Ohio (vector)
Custom layer I imported from a Shapefile containing lines (vector)

I am wanting to use these as an offline map in an iOS app.  Can I somehow combine these into a single mbtiles file?

I used a patch script I found on github, all tiles are there, but not the metadata for them (It only lists vector), I opened it in DBBrowser and can view the hillshade blobs. Is there a tool to combine these?


texts: Applying mask in ArcGIS Geostatistical IDW?

texts:
I'm finishing an analysis of approximately 850 archaeological assemblages in the Mediterranean, and I've successfully run IDW in Geostatistical toolbox (ArcGIS 10) to generate the results I like.  I'm using the Geostatistical toolbox as I'm interpolating the values of polygons, not points.  I'm trying to limit the extent of the results to the bounds of the Mediterranean Sea, however, but adding the mask to the appropriate dropdown box in IDW > Environments > Raster Analysis makes no difference to the extent of the interpolated results.
I tried Diffusion Interpolation with barriers, but that didn't provide the appropriate results visually, and I tried Extraction (Spatial Analyist > Extraction > Extract by Mask) as well.  The latter provided me with results within the geographic bounds I wanted, but appears to have reinterpreted my raster and is displaying the results in black-and-white - my original IDW results disappear.
At this point, I'm simply importing the resulting rasters into Photoshop and masking them there, but that can't be right.  Is there a way to mask IDW in the Geostatistical toolbox that I'm missing?


texts: How to use right-click to finish drawing a Polygon in OpenLayers?

texts:
I'm using OpenLayers 2.13  and the DrawFeature control. By Default, the Polygon is completed, when you double click.

Some of our clients are requesting that we allow the mouse right-click, to finish the polygon.

I tried adding a rightclick function to the OpenLayers.Handler.Polygon class, but it is never called.

I also confirmed that the right click event is being registered in the activate method of the OpenLayers.Handler class.

How do I override this behavior to allow the mouse right-click to finish the polygon?


texts:
I have a raster on ArcGIS and I want to add some contours. I use the ArcPy library to do that. This is the raster (part of a country) and the contours:

This is the raster with the contours:

and this is the same image, zoomed:

I am trying to do the same like the image below:

But the contours that are created are not in the shape I imagined. I just have a raster with different values and I want to add contours like in the last image.
Is there any tool that I can do that? In my script I use the simple arcpy.sa.Contour command. I just want to add contours of ground acceleration around an epicenter.


texts: ArcGIS Runtime SDK for .NET (Windows Phone) - ArcGISLocalTiledLayer Issue

texts:
I'm developing Windows Phone 8.1 client for our company ArcGIS maps. Online layers works fine, but I need to deliver offline maps (Map Tile Package - .tpk file) to users with bad connectivity. My problem is that not every user has device with Storage Card, and I try to save .tpk file to Internal Storage of Phone, i.e. Pictures Library folder.

I can succesfully read .tpk file from there to StorageFile, so the app has the permission to read, but I cannot initialize ArcGISLocalTiledLayer from this file. Even when I copy file to App Local Store, still cannot initialize ArcGISLocalTiledLayer from this file. I tried different path formats:

// Settings.LocalDataPath = "C:\\Data\\Users\\Public\\Pictures\\Camera Roll\\ArcGIS\\MS_PM_BaseMap.tpk"; - when loading from Pictures Library
// Settings.LocalDataPath = "ms-appdata:///local/MS_PM_BaseMap.tpk"; - when file is in Local Storage
// Settings.LocalDataPath = "isostore:/MS_PM_BaseMap.tpk"; - another try...

var tpk = new ArcGISLocalTiledLayer(Settings.LocalDataPath);
tpk.ID = "tpk";
await tpk.InitializeAsync();

on InitializeAsync() I get exception "Failed to open file. May not be a valid tile package.". When file is saved on StorageCard, initialization works fine. Is it possible to use .tpk file from any location in phone's internal storage?


texts: Image display/scaling in OpenLayers vs Google Maps

texts:
When I plot an image on OpenLayers, it looks like that image is stretched a little bit as compared to Google Maps. I was wondering if it is normal or there is something wrong with my code.

OpenLayers:

Google Maps (maps.google.com):

Following is the XML :

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
  <Folder>
    <name>Ground Overlays</name>
    <description>Examples of ground overlays</description>
    <GroundOverlay>
      <name>Large-scale overlay on terrain</name>
      <description>Overlay shows Mount Etna erupting
          on July 13th, 2001.</description>
      <Icon>
        <href>http://code.google.com/apis/kml/documentation/etna.jpg</href>
      </Icon>
      <LatLonBox>
        <north>37.91904192681665</north>
        <south>37.46543388598137</south>
        <east>15.35832653742206</east>
        <west>14.60128369746704</west>
        <rotation>59.1556640799496235</rotation>
      </LatLonBox>
    </GroundOverlay>

  </Folder>
</kml>

And following is my JavaScript code:

xml=jQuery.parseXML(oRequest.responseText)  //oRequest.responseText contains above xml

    groundOverlays=jQuery(xml).find('GroundOverlay');

    jQuery.each(groundOverlays,function(index, value){
        var groundOverlay=new Object();
        parentFolder=jQuery(groundOverlays).parent('Folder').find('name').text();

        groundOverlay.href=jQuery(value).find('href').text();
        groundOverlay.name=jQuery(value).find('name').text();
        groundOverlay.description=jQuery(value).find('description').text();
        groundOverlay.north=jQuery(value).find('north').text();
        groundOverlay.south=jQuery(value).find('south').text();
        groundOverlay.east=jQuery(value).find('east').text();
        groundOverlay.west=jQuery(value).find('west').text();
        groundOverlay.rotation=jQuery(value).find('rotation').text();

        var options = {
            'opacity': 1.0,
            'isBaseLayer': false,
            numZoomLevels : 20 ,
            //aspectRatio:50,

        };
        var image = new OpenLayers.Layer.Image(
            groundOverlay.name,
            groundOverlay.href,
            new OpenLayers.Bounds( groundOverlay.west, groundOverlay.south,  groundOverlay.east, groundOverlay.north),
            new OpenLayers.Size(1, 1),
            options
            );

        map.addLayers([image]);

    });


texts: How to use QGIS Time Manager sample data?

texts:
I've downloaded the Time Manager plugin for QGIS for Windows, version 1.7.1. I'm not sure whether I should really be asking a question on here, but since there doesn't seem to be an issue tracker on the new git, I want to post pictures, and I know the developer frequents this list, I figured I may as well.

My problem is I can't seem to make it work. Using the test data "tweets.shp", I can get an animation, however it is backwards, i.e. the screen is full at the start of the timeline, and tweets are removed at the time they are posted, leaving nothing at the end of the timeline. Graphics attached:

The settings I am using are:

start: t

end: none

index: tweets20120109101002200

timeformat: %Y-%m-%d %H:%M:%S

offset: 0

I installed it manually (due to firewall issues I think) into "\Program Files\Quantum GIS Wroclaw\apps\qgis\python\plugins\timemanager" - wonder if it's something to do with that.

Alternatively, is there a way I could debug the queries it is making to QGIS for each time point?


texts: Reduced/Increased size polygons with Buffering GDAL/OGR

texts:
In case somebody needs to do the buffering of geometries (using pure OGR).
I am trying to create reduced/extended polygons from an existing polygon using OGR in Python.
The abstract idea that I have is to use the "Buffering" function of OGR to achieve the goal that is dressed here and here (using GIS interfaces).
More precisely I want to create a buffer with negative/positive distance  and based on that create a new polygon with borders drawn inward/outward.
For this I used Create Buffer from GDAL/OGR cookbook as the template with a small addition of importing SRS from the input polygon to the output(buffer) polygon.
I have been suggested  to use SQLite dialect in order to achieve it. This means that there are certainly more options to get the same results.
Here is the code:
import os
import sys
from osgeo import ogr
from osgeo import osr

def createBuffer(inputfn, outputBufferfn, bufferDist):
    inputds = ogr.Open(inputfn)
    inputlyr = inputds.GetLayer()

    shpdriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outputBufferfn):
        shpdriver.DeleteDataSource(outputBufferfn)
    outputBufferds = shpdriver.CreateDataSource(outputBufferfn)
    # create the spatial reference system, WGS84
    #srs =  osr.SpatialReference()
    #srs.ImportFromEPSG(4326)
    srs = inputlyr.GetSpatialRef()
    bufferlyr = outputBufferds.CreateLayer(outputBufferfn,srs,
    geom_type=ogr.wkbPolygon)
    featureDefn = bufferlyr.GetLayerDefn()

    for feature in inputlyr:
        ingeom = feature.GetGeometryRef()
        geomBuffer = ingeom.Buffer(bufferDist)

        outFeature = ogr.Feature(featureDefn)
        outFeature.SetGeometry(geomBuffer)
        bufferlyr.CreateFeature(outFeature)
        outFeature = None

def main(inputfn, outputBufferfn, bufferDist):
    createBuffer(inputfn, outputBufferfn, bufferDist)

if __name__ == "__main__":
    inputfn = 'polygon.shp'
    outputBufferfn = 'polygonBuffer.shp'
    bufferDist = -1.0

    main(inputfn, outputBufferfn, bufferDist)


texts: Free GeoTIFF Overlay and Tiling (i.e. Indexing) Analyser?

texts:
Has anyone come across a shareware/freeware tool for analysing the overlays and tiling (i.e the internal indexing) of GeoTIFF images?

What would be excellent would be something like a GeoTIFF tuning advisor...

Back in my Intergraph days, we had such a tool.  It allowed us to have multi-gigabyte TIFFs being displayed in maps with subsecond performance.


texts:
I have a data set comprising a number of debris surveys around Australia. I wanted to aggregate these counts by grids of 10km by 10 km (approximately) and found a very useful post here: Aggregating points to grid using R

@whuber created a very clear way to graph the values, with the following code:

 count.max <- max(counts.m["Count",])
 colors = sapply(counts.m["Count",], function(n) hsv(sqrt(n/count.max), .7, .7, .5))
 plot(counts.m["X",] + 1/2, counts.m["Y",] + 1/2, cex=sqrt(counts.m["Count",]/100),
      pch = 19, col=colors,
      xlab="Longitude of cell center", ylab="Latitude of cell center",
      main="Event counts within one-degree grid cells")

Note that counts.m is a matrix with an X and Y value for the center points of the grids, and a statistic. In this case it is the total number of conflict events in Africa, but in my case it is the average amount of debris in a given grid.

I have re-created this code for my data set, but would now like to create a legend for the map. I've gotten as far as creating a fixed subset of boxes filled in with the appropriate colors, but can't get them to be circles or the appropriate size, as in the graph. Here is the code I have attempted thus far:

legend("topright", legend=c(1500,1000,500,25), pch=19,
   pt.cex=sqrt(c(1500,1000,500,25)/100),
   fill=sapply(c(1500,1000,500,25), function(n) hsv(sqrt(n/count.max), .7, .7, .5)))

Note that the original post (Aggregating points to grid using R) contains links to the data and code used to generate the count.max used above.


texts: ESRI Basemaps - World Imagery: Is it possible to bring certain imagery into foreground (change the viewing order)?

texts:
I have a basemap loaded in a map in which there are 4 overlapping images.  I would like to either turn off 3 of the 4 images or bring the one i want (marked with X n the screen capture) into the foreground.  Is this possible? Is there a workaround?

I tried un-checking lower resolution layers in the TOC but a messages pops up saying that thisis not possible.


texts: Storing geometry from a layer into a geopackage file with geotools plugin

texts:
How can I store a geometry to a geopackage file?

FeatureIterator<SimpleFeature> fi = getLayer().getFeatureIterator((Envelope) null);
GeoPackage geopkg = new GeoPackage(...);
geopkg.init();

final FeatureCollection<SimpleFeatureType, SimpleFeature> fc = new DefaultFeatureCollection(null, null);
 try {
    while (fi.hasNext()) {
        SimpleFeature sf = fi.next();
        Geometry geom = (Geometry) sf.getDefaultGeometry();
       ???
    }
} finally {
    fi.close();
}


texts: arcpy.da.How do I resolve arcpy.da.NumPyArrayToTable RuntimeError: An invalid SQL statement was used. [SELECT * FROM table]

texts:
I working on a python script to copy domains from one geodatabase to another (to rebuild replicas/staging environments).  I cannot figure out why I get the following error when it the script hits arcpy.da.NumPyArrayToTable(array, table)

Runtime error
Traceback (most recent call last):
  File "", line 70, in
  File "", line 65, in migrateDomains
RuntimeError: An invalid SQL statement was used. [SELECT * FROM table]

Below is my code

inputDomains = arcpy.da.ListDomains(inputDB)

tupLst = []

for domain in inputDomains:
    if domain.name == 'gFittingType':
        print('Domain name: {0}'.format(domain.name))
        print(domain.description)
        if domain.domainType == 'CodedValue':
            codeValues = domain.codedValues
            for val, desc in codeValues.items():
                tupLst.append((str(val),str(desc)))

        array = np.array(tupLst, dtype=[('code','S50'),('value','S50')])
        table = r'u:\scratch_workspace\scratch.gdb\table'
        arcpy.da.NumPyArrayToTable(array, table)
        arcpy.TableToDomain_management(table,'code','value',targetDB,domain.name,domain.description,'REPLACE')

Additional info (ArcGIS 10.3)


texts: How can I set specific zooms/resolutions in Openlayers 3 zoom control?

texts:
I want to have only specific zooms in my zoom cotrol so the map opens in zoom 3 , when click the zoom in it goes to 6 -> 12 -> 16 and back when zoom out.

How is this doable in OL 3 ?


texts: Are there trade offs between using a geohash vs using a quadkey as a spatial index?

texts:
QuadKey

Bing Maps uses a quadkey structure for their tiling scheme. Here is an overview of the concept http://msdn.microsoft.com/en-us/library/bb259689.aspx.

GeoHash

A geohash seems to be the more accepted representation in the open source word  http://en.wikipedia.org/wiki/Geohash

So what I am wondering is if there are any trade offs between the two when being used as a spatial index. They are both rooted in the quadtree concept and carry it's pros and cons, but are there benefits to using one over the other?


texts: Problem with QGIS Network Analyst (Basic)

texts:
I'm going through the QGIS training, Lesson: Network Analysis, and I've encountered a problem with Network Analyst. When I open the Shortest Path (Pt to Pt) tool, I can set the Vector layer representing network and Path type to calculate, but when I click on the "..." to set the Start or End Points, nothing happens.
I've searched on this site and Googled around and haven't found anything. I've tried it in both 3.10 and 3.8 and the problem exists in both.


texts: How to write Shapely geometries to shapefiles?

texts: Can someone demonstrate a simple way to write geometry data-structures from shapely into shapefiles? I am particularly interested in polygons with holes and linestrings. It would also be beneficial to stay away from arcpy (so osgeo, pyshp, etc. would all be better).


texts: Writing code with OFTInteger data type in OGR/GDAL?

texts:
I created the shapefile by using the feature field type: OFTInteger.

But when I open the shapefile in Arcgis (ArcCatalog) the data type is: Double (not integer).

Do you know how to fix this (I want it is Int) ?

Please see my codes:

OGRFieldDefn oSHToBDField( "SHBANDO", OFTInteger);
poSHPLayer->CreateField (&oSHToBDField);


texts:
I have to add raster data in geoserver using Web Coverage Service. I have raster data in geotiff format. The problem is that the geotiff image I have added to the geoserver is using WMS not the WCS facility. Please guide me how to add geotiff image using WCS service.

The 3 services of WCS are working properly but when I check the source code in javascript file, it shows the following code for the source of geotiff image

..../geoserver/webmap/wms?LAYERS=webmap%3Ayam2&STYLES=&FORMAT=image%2Fjpeg&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&SRS=EPSG%3A32643&BBOX=913482.6930837,6340905.95,746408.1930837,6421838.25&WIDTH=542&HEIGHT=330,

And in this code the image is using WMS service. Also, it uses GetMap request of WMS. I can't understand how it supports all 3 services of WCS perfectly even though uploaded geotiff image is using WMS service.


texts: PyQGIS - Start/Stop Editing with data Provider - Commiting changes

texts:
I want to edit a specific field for a QgsVectorLayer, within a QGIS application.

I am facing inconsistent results from my code, ie: changes are not alway commited. Please could someone point me to where it is wrong within this loop:

#vlayer is a QgsVectorLayer, imported from SQLITE DB
vlayer = QgsMapLayerRegistry.instance().mapLayer(layer)
caps = vlayer.dataProvider().capabilities()
polygonGeom = QgsGeometry( rubberBand.asGeometry() )
#if vlayer is capable of edition
if caps & QgsVectorDataProvider.ChangeAttributeValues:

    vlayer.selectAll()

    for feature in vlayer.selectedFeatures():

        if feature.geometry() == None:
            continue

        if feature.geometry().intersects(polygonGeom):
            fid = feature.id()
            idx = vlayer.fieldNameIndex('some field')
            attr ={idx:'some value'}
            print fid,idx,attr
            vlayer.dataProvider().changeAttributeValues({fid:attr})

    vlayer.invertSelection()

Here is a post of my console showing that some features are selected and should be edited, but again, some are not properly updated to the underlying database.

2 16 {16: 4}
12 16 {16: 4}
13 16 {16: 4}


texts: Using pyproj and Mercator: x_0=0 and y_0=0?

texts:
I'm having a hard time figuring out how to project correctly from a Mercator projection to LAT/Lon.

Basically, I have a netCDF file with the following variables:

LON0 = 42.0 => reference longitude for the conformal projection (°E)
LAT0 = -30.0 => reference latitude for the conformal projection (°N)
LONOR = 41.313 => longitude of point x=0,y=0 in the conformal projection (°E)
LATOR = -31.6997 => latitude of point x=0,y=0 in the conformal projection (°N)

I know that the point (476000,428000) is associated with geographical coordinates : ( 46.2559 , -27.8431 )

However, I can't figure out how to retrieve it by using pyproj transform object. Here is the code I'm using:

from pyproj import Proj,transform

projection_input = " ".join( ("+proj=merc +lon_0=42.0",
                              "+k=1 +x_0=0 +y_0=0",
                              "+ellps=WGS84 +datum=WGS84 +units=m",
                              "+no_defs")
                           )
inProj = Proj(projection_input)

projection_output = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
outProj = Proj(projection_output)

print transform(inProj, outProj, 476000, 428000)

Which gives me: ([46.27598102655299], [3.8677438926911805])

It has something to do with x_0,y_0 but I don't know what to put in it so
I can get the right transformation (projection). How can this be solved
by using the variables (only) listed at the top?


texts: Need to show a popup on selecting a layer option in leafletjs

texts:
I am trying to show a set of options using the layerGroup functionality to filter out the markers. I need to show a popup when an option is selected from a layer which will prompt the user to enter a number based on which i will be showing the markers. Is this possible as i am new to leafletjs?

In the below sketch, there' s a sample filter at the top right which i intend to show using the layers but the thing is on selecting the 'Location' option i have to show a popup that will prompt the user to enter a location number and on that basis the markers would should be populated.


texts: Finding if journey passes near certain address/location?

texts:
I am thinking to create a website where people can share car journeys. Users can register a car journey and other users can find a journey.

At the moment, when a user registers a journey, I save the starting and ending lat & lng of that journey in a database table called journeys. Other users can search a journey by entering their source and destination postcodes (postcodes a mapped to their lat and lng).

How can I know if a certain journey passes near the user searched destination postcode?

For example, a journey starts from location A and ends at location C but it passes location B. A user wants to find a journey from A to B.


texts: How to calculate x and y variance values?

texts:
I'm trying to define the best suitable cell size to use in the kernel density estimation in the Geospatial Modelling Environment (GME) for ArcGIS.

To calculate the best suitable cell size the following is advised (somewhere else):

  A cell size that is too small will result in a very large output
  raster ... that takes a long time to calculate. I suggest the
  following rule of thumb to calculate a reasonable bandwidth: take the
  square root of the x or y variance value (whichever is smaller) and
  divide by 5 or 10 ...

So, how should I proceed to the the x and y variance values?

All the best
HNL


texts:
I use osmconvert tool to merge maps. This tool merges 3 maps correctly, however if I merge 16 maps, then I see the following error:

  osmconvert Error: could not get 183500800 bytes of memory

This my command:

osmconvert.exe asia/armenia-latest.o5m asia/azerbaijan-latest.o5m
europe/belarus-latest.o5m asia/kazakhstan-latest.o5m asia/kyrgyzstan-latest.o5m
europe/moldova-latest.o5m russia-latest.o5m asia/tajikistan-latest.o5m
asia/uzbekistan-latest.o5m europe/bulgaria-latest.o5m europe/hungary-latest.o5m
europe/poland-latest.o5m europe/romania-latest.o5m europe/slovakia-latest.o5m
europe/ukraine-latest.o5m europe/czech-republic-latest.o5m -o=merged.o5m

Could you say what am I doing wrong? How can I solve this problem?


texts: Select All in the Python script tool doesn't work

texts:
When I click on the Select All, it doesn't work (using ArcGIS 10.3.1 Advanced). How can I fix this?

The Python Script for this model is listed below, from the online source rather than my creation.

# Import system modules
import arcpy

# Set the parameters
InputFeatureClass = arcpy.GetParameterAsText(0)
InputField = arcpy.GetParameterAsText(1)
InputValue = arcpy.GetParameterAsText(2)

arcpy.SetParameter(3, InputValue)

The parameter validation code is listed, also from the online source.

import arcpy
class ToolValidator(object):
  """Class for validating a tool's parameter values and controlling
  the behavior of the tool's dialog."""

  def __init__(self):
    """Setup arcpy and the list of tool parameters."""
    self.params = arcpy.GetParameterInfo()
    self.fcfield = (None, None)

  def initializeParameters(self):
    """Refine the properties of a tool's parameters.  This method is
    called when the tool is opened."""
    return

  def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parmater
    has been changed."""
    if self.params[0].value and self.params[1].value:
      fc, col = str(self.params[0].value), str(self.params[1].value)
      if self.fcfield != (fc, col):
        self.fcfield = (fc, col)
        self.params[2].filter.list = [str(val) for val in
                                        sorted(
                                          set(
                                            row.getValue(col)
                                            for row in arcpy.SearchCursor(fc,
                                                                          fields=col)
                                              )
                                           )
                                      ]
        if self.params[2].value not in self.params[2].filter.list:
          self.params[2].value = self.params[2].filter.list[0]

  def updateMessages(self):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""
    return


texts: Fill/remove holes in cloud shadow mask (Earth Engine)

texts:
I did a mask from image collection using SLC band of Sentinel. This mask represent the locations covered by cloud and shadow during the period picked.
Now a want fill/remove holes less than 1ha that showed up in mask. How I can do it?
I tried use focal_mean (), image.convolve () but these tools don't do exactly what I want. These tools increases the mask.

I tried use ee.Algorithms.FMask.fillMinima() but nothing happened (the layer stay equal the mask). I think this tool work to this situation. But I dont know how it works.
Follow my code:
// ID AREA AND TIME
var start = ee.Date('2020-05-15')
var finish = ee.Date('2020-05-30')
var tiles=[ '24LWL', '24LVJ']

//IMG COLLECTION with SCL (coleção S2A)
var imgCol = ee.ImageCollection("COPERNICUS/S2_SR")
  .filterDate(start, finish)
  .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 71)
  .filter(ee.Filter.inList('MGRS_TILE', tiles))
  .select(['SCL']);

//RECLASSIFIED
var reclassified = function(image) {
   return image.remap({
    from: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11],
    to: [0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1]
  });
};

var imgCol_reclass= imgCol.map(reclassified)
var imagem_2= imgCol_reclass.first()

//REDUCE
var img_reclass_reduce= (imgCol_reclass.reduce(ee.Reducer.product()))

//MASK
var mask = img_reclass_reduce.eq(1);
var maskedComposite = img_reclass_reduce.updateMask(mask);
var maskedComposite2= maskedComposite.visualize({palette: '#000000'})

//FILL MINIMA 1HA
var fillMinima= ee.Algorithms.FMask.fillMinima(maskedComposite, null, 100)

//FOCAL MEAN
var mean_scene = maskedComposite.focal_mean({
  radius:10,
  units:'meters',
  iterations: 1,
  kernelType:'square'
})

var focal_mean = maskedComposite.unmask(mean_scene)
var focal_mean= focal_mean.visualize({palette: '#DC143C'})

//ADD LAYER
Map.setCenter(-38.4872,-13.0156, 10);
Map.addLayer(focal_mean, {}, 'test_focal');
Map.addLayer(fillMinima, null, 'test_fillMinima');
Map.addLayer(maskedComposite2, {}, 'masked');


texts: Maplex Label Expression in Python Parser of ArcGIS for Desktop?

texts:
I am using a python code in ArcGIS to define the size of the labels of a point feature class. Therefore I am using a label expression (Maplex). The label size of the TEXTSTRING (string) should be defined by the TEXT_SIZE (double, with values like "3,5"). Generally the code is working, but not with special german characters (ü,ä,ö).

The code looks like that:

  def FindLabel ( [TEXTSTRING], [TEXT_SIZE] ):
      size = str([TEXT_SIZE])
      size_rec = str(float(size.replace(",","."))*4)
      size = size_rec.replace(".",",")
      return "<FNT size = '{0}'>{1}</FNT>".format(size, [TEXTSTRING])

Labels with these special characters are not displayed, labels with "normal" characters are displayed and sized correctly. Any ideas?


texts: Batch Saving Raster DTm and choosing type in QGIS 2.18

texts:
Follow on question from the following Previous question

batch saving 10m DTM layers using python console

I would now like to set the output of my Raster .tif files to Float32 as by default the .asc files I am saving are being converted to Float64 tif files.

Does anyone know how the previous quetions code could be modified to do this? Or if there is a similar way of doing this in QGIS3 I am open to using this software also as it is now in LTR with 3.4.5

Here is the Code being used-

myDir = 'C:\Users\david.laverty\Desktop\Save tif from asc/'
layers = iface.legendInterface().layers()
pipe = QgsRasterPipe()
for layer in layers:
   extent = layer.extent()
   width, height = layer.width(), layer.height()
   renderer = layer.renderer()
   provider=layer.dataProvider()
   crs = layer.crs().toWkt()
   pipe.set(provider.clone())
   file_writer = QgsRasterFileWriter(myDir + layer.name() + "_2017"+ ".tif")
   file_writer.writeRaster(pipe,
                           width,
                           height,
                           extent,
                           layer.crs())


texts: How can I build a custom address locator?

texts:
How can I geocode internal locations of, for example, an airport terminal via a custom address locator?

All I have is the footprint of the building itself.

I started by drawing out the footprints of the individual amenities within the building (which I'm able to do fine via the Editor tool). I also created an address file with store name, type, and concourse letter. Unfortunately, my building footprint does not reference any of these fields. I was hoping to geocode my "address" table onto the polygons I created.

Is there a reference to Esri lessons (for a custom address locator such as this)?


texts: Building FGDB support for GDAL 1.10 on Ubuntu

texts:
Because data.gc.ca is now distributing some of its data in FGDB format (example: Plant Hardiness Zones), I would like to be able to use it with gdal/QGIS. The compile/install is not going as planned, however.

These are the steps I've taken:

Enabled the source repository for ubuntugis-unstable.
Downloaded the Linux 64-bit File Geodatabase API 1.3, installed it in /usr/local/FileGDB_API, added the necessary ldconfig setting, and successfully built and run the test programs. (As suggested procedure here: [Qgis-user] How to compile gdal with ESRI FileGDB)
Installed the necessary source dependencies: sudo apt-get build-dep gdal
Pulled down the source: apt-get source gdal
Added the line --with-fgdb=/usr/local/FileGDB_API \ to the configure command in the debian/rules file.
From the source directory, issued dpkg-buildpackage -rfakeroot -uc -b (as advised here: How do I get the source code of packages installed through apt-get? - Ask Ubuntu)

The build goes well, but fails late in the package compilation stage, complaining about a missing Ruby 1.9 tree in the fakeroot filesystem:

…
mv /home/scruss/gdal-1.10.0/debian/tmp/usr/local/lib/site_ruby/1.9/* /home/scruss/gdal-1.10.0/debian/tmp/usr/lib/ruby/1.9/.
mv: cannot stat ‘/home/scruss/gdal-1.10.0/debian/tmp/usr/local/lib/site_ruby/1.9/*’: No such file or directory
make: *** [install] Error 1
dpkg-buildpackage: error: fakeroot debian/rules binary gave error exit status 2

As this is unrelated to FGDB, it looks like I've got everything else in place. But it also likely means that gdal wouldn't build anyway.

I'd rather stick with the ubuntugis packaging, as Ubuntu/Debian gets difficult with source packages built outside its control.


texts: How to create a new layer with specific data from a large data base of irrelevant data?

texts: I am a complete beginner to QGIS. I am trying to isolate an individual set of data from a very large data base to plot in another layer. I have isolated the piece of data by using the Selection by expression and entering the value I want. however, I am unable to duplicate that layer to just have the selected single data point I want in it. Does anyone know how to solve this basic issue for a newbie.


texts: Sorting string column in QGIS Print Composer?

texts:
Previously asked a the same generic question (Sorting string column Qgis), but would need to know a solution in Print Composer.
A added an atribute table in my print composer, and I need to sort the field.

I would like to order my field like "normal" but as this column is string and it doesn't work out. I got "Emplacement réservé 105";"Emplacement réservé 106"; .... before "Emplacement réservé 20".

What should I do? to have

Emplacement réservé 10
Emplacement réservé 25
Emplacement réservé 99
Emplacement réservé 105
Emplacement réservé 105


texts: Seeking Free Shapefile of European Countries?

texts:
I do not seem to be able to find a free shapefile of the European countries. This is the kind of map I am looking to populate with my own data using ESRI ArcGIS:

The Eurostat website only seems to have a shapefile of the region with country boundaries that also includes the NUTS divisions. Natural Earth shapefiles only seem to be available for the entire world, not for regions like Europe.


texts: Extracting {{bbox}} value of entire window in Overpass Turbo

texts:
I am using Overpass Turbo, quite a nice tool, but I want to extract the bbox values and use them in an osm-python-tools query.

Oddly I don't see any way to get the values of my current bbox, I don't see it in the UI anywhere, nor do I see it in the output of the Overpass QL query.

Maybe you can specify that it be output? Seems like it should be somewhere?


texts:
We are trying to enable Geodatabase (Oracle 18c) archiving (https://desktop.arcgis.com/en/arcmap/10.6/manage-data/geodatabases/geodatabase-archiving.htm)

Our architecture team recommends we separate out the tablespaces for primary tables and archival tables and use 'SAN (Block) - XIV' for primary and 'NAS (File)- NetApp' storage for the archive tables.

Now, since the data gets archived as users edit, will this have any performance impacts?

Also, archive data is what users will be utilizing to do 'time-lapse' animations and such using UI tools.

Will this storage impact the performance of such tools?


texts: Filtering area through defined time period using Google Earth Engine?

texts:
I tried the following code in the Google Earth Engine through the various countries which have different dry season. But I cannot understand how to iterate the code. can anyone generate iteration using this code over any country of the world?

//clip boundary for each and every collection
var clipBoundary = function(image)
                        {
                          return image.clip(countries.filterMetadata("ADM0_NAME","equals","India"))
                        }
//create a time band for linear regression
var createTimeBand = function(image) {
   // in the linear regression output.
  return image.addBands(image.metadata('system:time_start'));
};

 var startdate = ee.Date('2017-12-01');
  var enddate = ee.Date('2018-05-31');
var collection = ee.ImageCollection('MODIS/006/MOD13Q1').select('NDVI')
                        .filterDate(startdate,enddate)
                        .map(clipBoundary)
                        .map(createTimeBand)
                        .select(['system:time_start', 'NDVI'])
                        .reduce(ee.Reducer.linearFit())
                        .select('scale');
 print(collection)
  Map.centerObject(collection, 6)
 Map.addLayer(collection)


texts: Fill the color of polygons with some attribute value using ArcGIS

texts: I have a large set of polygon features where I need to add some attribute value in "Name"  field. Because it is a large set of polygons I want to have all polygons where I have add the value to be filled with some color and the empty ones to be uncolored. I am using ArcGIS 10.4. I have tried to create some new "Status"  field depending on "Name" field with value 0 for unpopulated polygons and value 1 for populated polygons. On that way I can set the symbology by "Status"  field. But I don't have idea how to do that.


texts: Tool/Algorithm for generalizing OSM-Files for printing maps

texts:
Is there a tool or algorithm to (topological) generalizes OSM-shapes? I would like to build a printable map in different scales.

In the image below you can see the same divided road (yellow) at different scales.


texts:
I'm looking to create XYZ tiles for use in leafletjs with nodejs. Previously I've been creating my tiles using this python library (gdal2tiles).

I was wondering if anyone had yet written a nodejs equivalent?

I'm looking to take a single png image and specify the zoom range for the tiles to generated at, min zoom probably defined by log2(max(width, height)/tilesize)

Ideally each zoom range could be created by an independent process as I'd like to spin up a separate AWS lambda instance for each zoom level and upload the tiles to an S3 bucket in parallel.

How do I generate these tiles in JS rather than Python where I lack the knowledge to edit?


texts: Reading out which map theme is used in QGIS 3

texts:
I'm using QGIS 3.2.2 and Python. I would like to read out which map theme I currently use.

I know the following method provides me a list of defined map themes ('Theme1', 'Theme2'):

QgsProject.instance().mapThemeCollection().mapThemes()

I don't know any method which provides me with the map theme name used itself.

I tried to use iface.mapCanvas().theme() but it always returns an empty string.


texts: Reclassifying Vector Field using python parser in field calculator

texts:
I am attempting to reclassify a field in ArcMap using python within field calculator.

Below is screenshot of how it is set up. when I run the reclassify I get an error, and my prepossessing results give me an error in syntax. Also attached is the geprocessing results. I am sure this is a very simple fix, I am just very new to python.


texts:
I have some spatial processing in R that I want to do.  After selecting out the appropriate polygon (clip1) I want to crop the raster on that polygon with the command clip1 <- crop(imported_raster, extent(single))

However, occasionally I get an error Error in .local(x, y, ...) : extents do not overlap because clip1 is not within imported_raster

> extent(imported_raster)
class       : Extent
xmin        : 143229
xmax        : 243229
ymin        : -564760.8
ymax        : -464760.8
> extent(single)
class       : Extent
xmin        : -133957.8
xmax        : -133157.8
ymin        : 165319.3
ymax        : 166119.2

Is there a simple/efficient way of tryCatch'ing this? My first attempt (below) is obviously not working.

> try({
+      clip1 <- crop(imported_raster, extent(single))
+      clip2 <- rasterize(single, clip1, mask=TRUE)})
Error in .local(x, y, ...) : extents do not overlap
>
>  tryCatch({
+      clip1 <- crop(imported_raster, extent(single))
+      clip2 <- rasterize(single, clip1, mask=TRUE)})
Error in .local(x, y, ...) : extents do not overlap


texts: Disabling OGR_GPKG_FOREIGN_KEY_CHECK option in QGIS

texts:
I am trying to load a bunch of vectors and rasters in a geopackage file.
I was able to load few vectors but then all of a sudden the following warning popped up preventing me to load any more vectors:

Failed to import some vector layers!
Opening of data source in update mode failed (OGR: pragma foreign_key_check on 'F:/link/to/my/geopackage.gpkg' failed. You can disable this check by setting the OGR_GPKG_FOREIGN_KEY_CHECK configuration option to NO)

I couldn't find a way to set that configuration or to manually override it.

How can I overcome this issue that prevents me from loading any more vector files?


texts: Publishing maps and services to ArcGIS Server 9.3 when created by ArcGIS Desktop 10.0?

texts:
We have license of  (ArcGIS desktop10) for our users.

We use it for create maps, when we want to publish that maps and services by licensed ArcGIS for server 9.3 , we first  convert MXD  (version 10) to  version (9.3) then publish it, The problem is we miss many things (Resolution of layers[Satellite image, features(line,...],some tools like (pan, identify tool,...) not work as well, and  some data not display!!
(NOTE: we don't have ArcGIS for server10)

My questions are:

What should we do to get a maps and services with more accuracy and clarity ?
Can we do any thing to  publish maps and services (that created by ArcGIS Desktop10) by ArcGIS for Server 9.3 directly without converting versions?


texts: Tool in ArcGIS for stream flow analysis?

texts: I'm just wondering if there is any tool in ArcGIS that can be used to analyse stream flow measurements. For example, estimating runoff, baseflow and baseflow index based on stream flow measurements.


texts:
I would like to compare the data of two geodatabases, in SDE or a file geodatabase,  to detect the changes occured between them , i've tried to compare the Xml workspace document of both databases using notepad++ "compare plugin" , but it's not easy to check and interpret the differences.

i'm looking for a better solution.


texts: ArcGIS 9.3.1 VBA: Unable to loop through datasets using 'ILayerGeneralProperties' when raster catalog is present

texts:
I have a simple project which involves the use of several buttons in order to make various datasets visible depending on which category is required at the time.

It works fine using the popular GetLayerByName function (http://bit.ly/Ha2iGK) but unfortunately there are duplicate layer names scattered across the mxd.  To get around this I modified it slightly to look at both the layer name and the layer description, but as raster catalogs are unable to store a description (this is unfortunately not in a geodatabase in case that makes any difference) the code falls over.

I've got a separate routine that clears all of the layers once the button is clicked and then it turns on the layers listed under another sub using [layername], [layerdescription]:

    Set pLayer = GetLayerByName("Wet (483)", "L6_Dry")
If Not pLayer Is Nothing Then
    pLayer.Visible = True
End If

Does anyone have any ideas or suggestions on how I can get it to skip/ignore the catalogs?  I've never had much dealing with them inside VBA so any help at all would be greatly appreciated.

Public Function GetLayerByName(sLayerName As String, sLayerDesc As String) As ILayer
Dim pMxDocument As IMxDocument
Dim pCompositeLayer As ICompositeLayer
Dim pMap As IMap
Dim l As Long
Dim m As Long

Set pMxDocument = ThisDocument
Set pMap = pMxDocument.FocusMap

For l = 0 To pMap.LayerCount - 1
    Set GetLayerByName = LayerByName(pMap.Layer(l), sLayerName, sLayerDesc)
    If Not GetLayerByName Is Nothing Then Exit Function
    Next l
End Function

Private Function LayerByName(pLayer As ILayer, sName As String, sLayerDesc As String) As ILayer
Dim pReturnLayer As ILayer
Dim pCompositeLayer As ICompositeLayer

Dim l As Long

Dim pLayerGenProp As ILayerGeneralProperties

'#### This is where the error appears!####

Set pLayerGenProp = pLayer

Dim thestring As String
thestring = pLayerGenProp.LayerDescription

If UCase$(pLayer.Name) = UCase$(sName) And UCase$(thestring) = UCase$(sLayerDesc) Then
    Set LayerByName = pLayer
    Exit Function
End If

If Not TypeOf pLayer Is IGroupLayer Then Exit Function

Set pCompositeLayer = pLayer
For l = 0 To pCompositeLayer.Count - 1

    Set pReturnLayer = LayerByName(pCompositeLayer.Layer(l), sName, sLayerDesc)

    If Not pReturnLayer Is Nothing Then
        Set LayerByName = pReturnLayer
        Exit Function
    End If
    Next l

End Function


texts: I'm using QGIS 3.4 and I want to know how to convert a column datatype from string to datetime in python. I found that we can do that by to_datetime from the calculator but I want to do it from Python. Here is my data like:


texts: How to convert a track of points into walking directions?

texts:
Imagine you have a list of points (latitude + longitude) you have previously recorded with a GPS receiver. Now, you need to convert it to a textual representation of the route you followed. Example:

Start at Redmond Avenue 12
Walk 200 m
Turn left to Another St.
Walk 100 m
Stop at Another St. 10
End

It doesn't really need to be very detailed.
The first idea that comes to my mind is to use an external service (google, yahoo...) to reverse geocode each point. Then apply a simple algorithm to write the route when I detect a change in the street name.

Any tip about what is the best algorithm to use or any existing (open) implementation is welcome.  It does not seem to be trivial.


texts: Using ArcPy to join feature class to Oracle Spatial View?

texts:
I am trying to join an Oracle table to a feature class but I keep getting the error "000339 : Input  does not have OIDs."

The view has a unique incrementing field already, but I am unsure how to specify with Python which field is the key?

Here is my code:

# Import arcpy module
import arcpy
from arcpy import env

#Set Variables
parent = "Database Connections\\db.sde\\db.DBO.layer"
child = "Database Connections\\db.sde\\db.DBO.layer2"
codeList = ['SC21', 'SC22', 'SC23', 'SC24', 'SC25', 'SC26', '2C27', 'SC28', 'SC29', 'SC30', 'SC31', 'SC33', 'SC34', 'SC35', 'SC36', 'SC37']
parentView = "layer"
tableView = "connectionFees"
table = "Database Connections/oracle.sde/oracleDB.oracleTable"

if arcpy.Exists(child):
    arcpy.Delete_management(child, "FeatureClass")

#Make a feature layer in order to perform a selection
arcpy.MakeFeatureLayer_management(parent, parentView)

#Select all features that match the codes from the code list
for code in codeList:
    arcpy.SelectLayerByAttribute_management(parentView, "ADD_TO_SELECTION", "Code = '"+ code+"'")

#Copy the selected features
arcpy.CopyFeatures_management(parentView, child)

#Make table view of Oracle table since it does not contain an incrementing ID field
arcpy.MakeTableView_management(table, tableView)

#Join the gdb layer to the Connection Fees Table from oracle DB
arcpy.JoinField_management(child, "field1", tableView, "field2")

UPDATE:

I have tried to swap the line

arcpy.MakeTableView_management(table, tableView)

For

arcpy.MakeQueryLayer_management("Database Connections/db.sde", tableName, "SELECT   \n         field1,\n         field2,\n         field3,\n         SUM (amount) total\n    FROM layer\nGROUP BY field1\n", "field1", "", "", "")

While this is successful, when I run

arcpy.JoinField_management(child, "field1", tableName, "field2")

I get a generic 999999 error.. Any thoughts?

UPDATE #2:

Turns out there is a bug in ArcMap. When I make the Query Layer, ArcMap automatically adds a percent sign to the beginning of the layer name. This is preventing the Join, causing the 999999 error. My best workaround is to copy the features into the goedatabase and then perform the join. I will post a more detailed answer below.


texts: Cannot insert data into my feature class

texts:
I'm using Arcgis 10.1 and I created my script for inserting data in my feature class but I got this error message:

Traceback (most recent call last):
  File "C:\SIG_AS\Shapefiles\Projet2\MyScript\Insert_chantiers_gp.py", line 35, in <module>
    with arcpy.da.InsertCursor( inGeodatabase, fieldsToUpdate) as cursor:
RuntimeError: cannot open 'Chantiers.shp'

And this is my script for inserting my data into my feature class:

# Ajout de champs
import arcgisscripting
import arcpy
# Définir l'environnement
arcpy.env.workspace = "C:/SIG_AS/DB/Chantiers.gdb"

# creation de l'objet gp (arcgisscripting est l'ancêtre d'arcpy et devait être créé à chaque fois)
gp = arcgisscripting.create(9.3)
# Paramètres d'entrée
inDescription = gp.GetParameterAsText(0)
inRestriction = gp.GetParameterAsText(1)
inNumero = gp.GetParameterAsText(2)
inNomVoie = gp.GetParameterAsText(3)
inDate1 = gp.GetParameterAsText(4)
inDate2 = gp.GetParameterAsText(5)
inLat = gp.GetParameterAsText(6)
inLong = gp.GetParameterAsText(7)

# Variables de champs
inGeodatabase = "Chantiers.shp"
ObjectidField = "OBJECTID"
ShapeField = "SHAPE"
descriptionField = "Desc_travaux"
restrictionField = "Restriction_trafic"
numeroField = "Numero_civique"
nomField = "Nom_voie"
date1Field = "Date_debut"
date2Field = "Date_fin"
LatField = "Latitude"
LongField = "Longitude"

#Tuple pour mise à jour des champs
fieldsToUpdate = ( ObjectidField, ShapeField, descriptionField, restrictionField, numeroField, nomField, date1Field, date2Field, LatField, LongField)
# Curseur d'insertion
with arcpy.da.InsertCursor( inGeodatabase, fieldsToUpdate) as cursor:
# Insertion de ligne aux attributs considérés
    cursor.insertRow((inDescription, inRestriction, inNumero, inNomVoie, inDate1, inDate2, inLat, inLong))


texts: Standard Method to calculate Average Moving Speed from Times & Coords?

texts:
...I wasn't sure if this question was best suited to GIS, SO or Math, so after reading GIS Topics I decided to ask here, but I'll relocate if necessary...  :-)

I'm working on automating various analyses of movement and activities based on various sizable (50K→500K) sets of datetime/latitude/longitude data.

I'd like to determine "Average Moving Speed".

Parts of the recordset will have little or no movement
Other parts will have movement of varying speeds (ie., walking or driving), and,
A few records will (frustratingly) contain what I assume are GPS errors: Sudden "hops" to an unrealistic distance, returning to the previous location in the following record.

Is there a standard method of determining Average Moving Speed?

Perhaps a minimum distance that is considered "actual movement" as opposed to GPS resolution or reporting issues, or
Maybe there is a standard formula or algorithm that I haven't found yet.

I'm comfortable with coding and I've brainstormed and experimented with a few methods but each seems to have it's own issues.

Possible methods:

Skipping every n records  (Surely the worst option.)
Skipping records between "whole minutes" or "whole 3 minutes" etc (Some legitimate movement would be missed.)
An interpretation of Rolling Average.

I figure that once I can extrapolate Average Moving Speed, I will be able to deduce Activity Type -- ie., whether the movement was most likely walking, driving, flying, idle, etc... but I'd be interested to hear thoughts about that as well.

From what I understand Rolling Average would be, instead of looking at the current record, take an average of n records before and n records after the current record.

If that's the case, there could still be various interpretations:

Take separate averages of the n + n + 1 Lon's & Lat's and use those numbers to calculate movement
Averages of each n + n + 1 records elapsed time, divided by Averages of each n + n + 1 records' distance traveled
I will try out any ideas you may have and report my findings back here.

An example of the data (with distance [meters] and time [seconds] added, at individual-record level):   (and this link to a downloadable 200-record CSV sample)

Some places I've looked for insight (without success, unless I missed something):

How to best measure average speed
How can I calculate the average speed in various legs of a bus route
Delayed Moving Average
Wikipedia: Moving Average


texts: Why does raster's cell size change due to rotate and georeference using Georeference toolbar?

texts: I just found out that the cell size of raster change after I rotate or georeference it with Georeference. Is there any explanation why it works that way?


texts:
I have three layers:

building data (polygon) attributed with parking lot capacity
on the road parking data (line), also contains parking space data
car parks (polygon) includes parking space capacity

I would like to create a raster, the raster must be 50 x 50m resolution, in every 50m square I want to know the parking capacity in it, how do I divide my city into pixels and calculate this parking space capacity in the pixels?  I am using ArcGIS 10.1


texts:
I used the Geographically Weighted Regression (GWR) to find out which independent variables lead to crime and how they vary over the space.

After a literature review, I ran a global Ordinary Least Squares (OLS) model to find the most suitable set of independent variables, then used the local GWR method. For both the OLS and GWR model I then conducted the Moran's I for the standardized residuals, to be sure that they are not spatially autocorrelated. I used ArcGIS for the analysis.

When I investigated the coefficient values within the different regions in the study area, I came to an issue that is unclear for me:
For example I have included the risk factor bus stop, given as density per square kilometers, to explain burglaries (given as crime rate per 100 000 inhabitants).

The values for C1_BusDen (Coefficient 1 for bus stop density) ranges from 21 to 39.
One region  has a very low value with a crime rate of 1100 and bus stop density of 4.5.
Another region has a very high value with a crime rate of 1250 and bus stop density of 5.07.

If the density of bus stops is high in both regions and there is also a similar crime rate, shouldn't be the relationship  (value) also the same - or at least similar, but not in the lowest and highest class?
Also, the region with the highest density of bus stops has at the same time the highest crime rate - the GWR coefficient value is only 26 - shouldn't it be the highest?

Do I interpret the result in the correct way, because it would make sense that a higher value in the independent variable and a higher value in the dependent varialbe lead to a stronger relationship.


texts: Creating Custom Shapefile that can be snapped in location with no dragging using QGIS?

texts:
I am trying to create a custom shapefile consisted of 16 waypoints orientated around an origin point. the 16 waypoints are to be 2m apart from the central waypoint (like an axis, with 0 being the central waypoint and 4 waypoints in each N,S,E,W direction). I have used geopackage to create this however i need to ensure they are exactly 2m apart. Is there a way to ensure this.

I also want to be able to use this shapefile in different locations, however do not want to drag the shapefile over long distances.

Is there a way to create this shapefile?

I can elaborate if needed - imagine placing an axis over your house and every 2m represents a sample point. Then you want to place the same axis over your friends house which is ~1,000km away.

Can we zoom to the friends house and snap the axis on easily?


texts: How to sum pixel values across featureCollection polygons with the same property?

texts:
I'd like to use the hansen forest cover dataset and the Resolve ecoregions database to sum forest areas across biomes.  The Resolve dataset has many polygons, but the "BIOME_NAME" column is repeated across many of them, such that there are only 13 biomes across the globe.  I would like to sum the forest area in each biome.

I'm sure there's some way to do this with reduceRegions, but I'm not sure how to confine the analysis across biomes, and not for each and every polygon.  At present, the code below results in an error when calculating the final dictionary.

https://code.earthengine.google.com/2a60e1ddacbafdcb70b1bffb59811812

var gfc2018 = ee.Image("UMD/hansen/global_forest_change_2018_v1_6"),
resolve = ee.FeatureCollection("RESOLVE/ECOREGIONS/2017");

var fc2000_pct = gfc2018.select(0).divide(100)

// Get scale (in meters) information from band 0.
var b0scale = gfc2018.select(0).projection().nominalScale();
print('Band 0 scale: ', b0scale); // ee.Number

var forest_area_per_pixel = fc2000_pct.multiply(b0scale)

// now to aggregate by biome...
var forest_area_per_biome = forest_area_per_pixel.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: resolve.geometry(),
  scale: b0scale,
  maxPixels: 1e9
});

print(forest_area_per_biome);


texts: ModelBuilder Iterations and Data Storage?

texts:
I'm new to ModelBuilder.

I have a large series of rasters divided into two categories, let's call it Type A & Type B, a pair for each of a few hundred date/time snapshots. The date/time for the snapshot is in the file name (roughly in the middle with a standard structure), and the Type is in the filename (near but not at the end). They're all in a folder together. I need to calculate a new raster for each date/time snapshot using the Type A raster & Type B raster to do so.

I know how to iterate rasters, and I know how to use raster calculator but how can I can store a raster while I'm iterating so that I can use both A & B?

The way the file names are, Iterate Rasters will run through Type A then Type B for each date before moving on to A for the next date, which I think is exactly what I want.

I just don't know much about how the variables work in ModelBuilder or ArcPy or even Python for that matter, so this is where I'm stuck.

I figure there's a few different ways to do it, and two that I've tried and failed (due to my limited knowledge of how to set it up) are:

Storing Raster A in a variable and when Iterate Rasters hits raster B, having it realize it's B (through a count variable or through a substring in the file name, etc.) and using the pulled-up raster B and stored raster A to calculate a new raster, then clearing both my raster storage and whatever other variables I have.
Using a built-in function like Collect Values and clearing that every two iterations. I'd have to know how to clear that and also how to calculate a raster using a reference to its position or whatever in Collect Values.

Can anyone explain in more detail (variable type, path, setup, code for things, etc.) how I could do either of those?

Can anyone tell me the code necessary in a Calculate Value function (or whatever else) to extract the characters between two known strings in a filename, that would give me the date?


texts: Reading cell values of Aspect map using ArcPy?

texts:
I have generated an aspect map from a DEM using ArcGIS and I want to create a CSV file containing the cell values using the arcpy library. I have done reading raster maps using the arcpy library previously with "RasterToNumPyArray" function. But aspect map cannot be read like that.

How to read the cell values in aspect map and generate a CSV file?


texts:
I am clustering rainfall data in my country using QGIS. Unfortunately, there are errors that I encountered during the process of interpolation. I already have a vector file with multiple attributes (Jan rainfall - Dec rainfall).
In the K-mean clustering method, I disable the "Normalize Attributes" option since I am afraid that it will affect the result of my project. But an error pops up saying "Problems during clustering... type other than float or double not supported". But if I enable it, the clustering proceeds.

In the Hierarchical-2 method (unsupervised). I do not know what to put on the following boxes since every time I run the plugin, it creates thousands of clusters.


texts:
Is the "sdetable" command available with ArcGIS Workgroup, or only the Enterprise edition?

We have ArcGIS Server Workgroup 9.3.1 and want to update some of the feature class tables with data from another system using SQL UPDATE commands.  From taking a look at the documentation, we must first create a multi-versioned view, and then run the SQL UPDATE commands against this view.  The problem is that we do not have the "sdetable" command required to create the multi-versioned view.  We think that these command line tools may be available only in the Enterprise edition.

Is there another way that we can perform basic data updates with SQL commands in the Workgroup Edition?  I understand that the correct way to do this is not by editing the default version directly.

EDIT:
Our data is versioned for different uses, including some field applications.  We were planning on creating a "sql-updates" version just for making updates to attributes in a feature class.  If I use SQL UPDATE commands to update this version, my concern is that the changes will not get replicated to the other versions, or the data could become corrupt.  This seems like it would be bypassing delta table entries.

From what I have read, we should use "sdetable" to create a multiversioned view, and then use stored procedures create_version, set_current_version, edit_version (1 - start session),  run my updates, edit_version (2 - close session), and finally delete_version (after reconciled).  I see these stored procedures in the database, but can't do the first step of creating the multiversioned view.  Maybe it would be possible to perform these steps on the "sql-updates" version without creating the view.

EDIT 2:
Since the original question has been answered, I opened a new question:
How can we use SQL to update feature class attributes in a versioned geodatabase in ArcGIS Workgroup


texts: ArcGIS Routing on Android with Network Analyst module

texts:
I am currently developing an Android App for navigation purposes. It is based on the ArcGIS Android API. I have an ArcGIS server (with network analyst module) with all the needed data running.

At the moment the App displays a map with the current position (point A) and has a search form, where you can search for a building (point B).

The next step would be to display a route from "point A" to "point B". Does anybody know if there is a way to do that natively with the Android API?

Also it would be possible to compute the routing information on a Java web service and then send it to the app, which only has to display it.

I cannot find any sample code on the ArcGIS resource page.


texts: Classified geologic maps in QGIS (3.12.2)

texts:
I'm really new to QGIS and trying to understand how to import and classify a geologic map. For example, here is a map from the AZ geological survey: http://repository.azgs.az.gov/uri_gin/azgs/dlio/1737

I've downloaded the GDB for the map and I can add the polygons etc, but they are all imported as a single category. Is there some way to get QGIS to create the symbology to look like the finished PDF? Please don't tell me that I need to pick each color by hand.


texts: Splitting shapefile into many shapefiles with open source?

texts:
I want to create a shapefile for each of the 72,000 census tracts in the United States that includes the census blocks in that tract using open source software.

I will probably start with a state shapefile that includes all of the blocks for that state and merge it with a GDB file. This will give me a list of blocks that are in the census tract.

But then I'm not sure how to split the state shapefile into smaller ones.

How can I do this?

Is there a way to automate this using open source software?

If so, what tools should I use?

I'm learning QGIS.


texts:
i'm new using QGIS (my version is 2.14.9) and received multi ortophotos (Supposedly Georeferenced in Datum 73) and so i'm trying to put them all  together using "Build Virtual Raster" following the Procedures described in Joining several raster files using QGIS?
but i get the following error menssage:

  "Warning 6 gdalbuildvrt does not support ungeoreferenced image"

What i did was "add layer" for which one of the orthos and sellected Datum 73 and then, when using the build virtual raster, instead of using "Choose input directory instead of files" i had choosen the option to use the layers i had added.

From what i can understand QGIS is not assuming that the ortophotos i have are georreferenced, but i don't know what i'm missing, can anyone please help?

PS: i have checked if the ortophotoss are georrefereced using AUTOCAD CIVIL 3D and, when choosing Datum 73 everything was fine.


texts: arcpy.env.overwriteOutput = True appears not to be working?

texts:
AOI_test = "AOI_test.shp"
arcpy.env.overwriteOutput = True

for i in loc_list:
    AOI = "AOI_" + i + ".shp"
    arcpy.Select_analysis("locality", AOI, " \"LOCALITY\" = '%s'" % (i))
    arcpy.Merge_management([AOI_test, AOI], AOI_test)

This code isn't working for some reason. Have I understood this code incorrectly, or should I be able to overwrite the merge output FC (AOI_test)  repeatedly, now that I have specified that this is acceptable?


texts: Finding neighbors for non adjacent polygons in R

texts:
I'm looking for a method or function in r-project that create a neighborhood matrix from a group of polygons that don't touch each other, like polygons of set of blocks separated by streets in a town.

I don't have a sample shapefile.


texts: Selecting points to create line in ArcMap?

texts:
Is there way in ArcMap 10.x to select a bunch of points and create a line between all the points?

There is a tool called Connect the Dots for MapInfo Pro that allows you to select any number of points and it will create a line connecting all the points.


texts: Is there a tool that shows dimensions during editing?

texts:
Is there a plugin or tool providing following functionality:

I select 2 points (or vertices) - it shows the length between this points.
When I edit the point (or move vertex) the length automatically updates (allowing me to easily move point to the distance I want).

I know there is dimensioning plugin but I think it is not what I want and it is not working well in 2.6.1.


texts: 1 - 1 Composite relationship class set to notify "both" not updating origin when editing destination

texts: I have two spatial feature classes and want to have the spatial edits reflected in both feature classes when editing either the origin or destination. I have set up a composite relationship with a 1-1 relationship and notifications set to both. When editing the origin the destination changes with it ( I believe this is default with composites set to forwarding). Problem is, when editing the destination, the edits are not reflected in the origin. I thought if set to both this would work out.


texts: ArcGIS Network Analyst tools are not using historical traffic data to calculate travel time

texts:
I have successfully created and built my network dataset with Historical HERE traffic data. However, Network Analyst tools are not using my traffic data to calculate travel time. No matter what time of the day I chose to run the tools, it only came back with the "free-flow" (or the Fallback attribute) travel time.
Here are some screenshots of my ND settings:


texts: How to extract slope value for a geometry point in GEE

texts:
I want to extract slope values for various points in a certain region. Below is my current code which visualizes my region of interest on the map viewer

// dataset
var srtm = ee.Image('USGS/SRTMGL1_003');

// slope
var slope = ee.Terrain.slope(srtm);

// region of interest represented as a Geometry Rectangle
var roi = ee.Geometry.Rectangle([35.86, -16.64, 34.26, -15.21]);

// clip image dataset to roi
slope = slope.clip(roi)

var visualization_params = {min: 0, max: 60, palette:'white,red'};
Map.centerObject(roi);
Map.addLayer(slope, visualization_params, 'slope');

My goal is to derive the slope values for different coordinates within the rectangle, whereby i can map a point (e.g. 34.26, -15.91) to its slope value, then export the mappings.

I don't know how to go about it, and would need some help.


texts:
What is the best definition of Geographical Information Systems?

Has the definition changed at all since I was taught it 10yrs back?
The technology has changed a lot and it is more powerful than it used to be, but has that changed the overall definition of what it is?

e.g. Is stuff like Augmented Reality a separate discipline altogether, or does it also merge into the world of GIS?


texts: Using ConvertToMSD in ArcPy to report an error

texts:
I have two questions:
1. Is the mxd file generated by arcgis10.4 not converted to msd file by this method?
2. If can, my code is like this

import arcpy
mxdPath = "F:/mm.mxd"
msdPath = "F:/test.msd"

mxd = arcpy.mapping.MapDocument(mxdPath)
print "---------------------------mxd open success-----------------------------------"

df = ""
frames = arcpy.mapping.ListDataFrames(mxd, "图层")
if len(frames) == 0:
    frames = arcpy.mapping.ListDataFrames(mxd, "Layers")
df = frames[0]
print "---------------------------------mxd to msd ing-----------------------------"
arcpy.mapping.ConvertToMSD(mxdPath, msdPath, df, "NORMAL", "NORMAL")
print "---------------------------------mxd to msd success-----------------------------"

the error is:

  AttributeError: 'str' object has no attribute '_arc_object'

all is this

File "f:/code/pythonTest/demo.py", line 17, in <module>
    arcpy.mapping.ConvertToMSD(mxdPath, msdPath, df, "NORMAL", "NORMAL")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\mapping.py", line 189, in ConvertToMSD
    result = mixins.MapDocumentMixin(map_document).convertToMSD(out_msd, data_frame, msd_anti_aliasing, msd_text_anti_aliasing)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\arcobjects\mixins.py", line 895, in convertToMSD    map = self._mxd._arc_object
AttributeError: 'str' object has no attribute '_arc_object'


texts: How to determine projection parameters when customizing a projection

texts:
I am trying to customize an Albers and a Hotine Oblique Mercator (HOM) projection to minimuze distortion in the region I am analyzing. The region extends from about 51 to 62 degrees latitude, covering an area about the size of the Ukraine. The region is oriented NW - SE.

I want to make sure I am using correct methods for determining the two projection parameters: lat/long of the projection center, and centerline azimuth.  I am using ArcMap v10.  Here's the procedure I've followed so far:

Created a single polygon that defines the analysis region (by, generally,
creating a convex hull around the extent of the watersheds covering
the region).  This polygon is the area I am customizing the
projection for.
Projected the polygon to Geographic/NAD 83.
Used Jeff Jenness' Tools for Graphics and Shapes
(http://www.jennessent.com/arcgis/shapes_graphics.htm) to determine
the polygon's center of mass on the GRS80 spheroid.  The resulting
coordinates are what I used for the "projection center" parameter.
To determine the centerline azimuth, I first projected the polygon
to an azimuthal equidistant projection, specifying the projection
center at the coordinates determined in Step 3.
Then I drew a polyline (in the azimuthal equidistant projection),
snapped to the projection center point, representing the directional
trend of the region polygon.  To get the azimuth at the projection's
center, I used Jeff Jenness' Tools for Graphics and Shapes to
determine the beginning azimuth of the geodesic curve at the central
point.
For the Albers projection I am using the longitude for the
projection center, as determined in step 3. Am also using the awesome spreadsheet
created by Bill Huber (http://forums.esri.com/Attachments/34278.xls)
to determine where to place the standard parallels to minimize the
scale distortion within the polygon region.

If needed to know, I am using ArcMap's version of the HOM that uses a central line defined by one point at the projection's center and its angle of azimuth. ESRI calls this the "Hotine_Oblique_Mercator_Azimuth_Center". In EPSG, I believe this is the Oblique
Mercator, Hotine Variant B, EPSG method code 9815.

I am hoping there are some projection experts out there that can tell me if the above procedure, especially Steps 3 and 4, is a correct way to determine the needed projection parameters. Am I on the right track? Is it correct to be determining the center on the spheroid and the angle of the geodesic from the center point (instead of a "2d" geometric center and azimuth)?

I hope the problem description was clear. I'm eagerly looking forward to any answers, tips, discussion, etc.!


texts: My NDVI values are way out of range - GEE

texts:
I am supposed to collect NDVI values for some states in the USA. I understand NDVI is supposed to be between -1 and 1, but my values are in the 1000s. Can someone please help me with the interpretation of the values I get or how to get it right? Here is the code I used:

var startDate = '2000-01-01';
var endDate = '2019-05-01';

var collection = ee.ImageCollection('MODIS/006/MOD13Q1').filterDate(startDate,endDate).select('NDVI');

var reducers = ee.Reducer.min().combine({
              reducer2: ee.Reducer.max(), sharedInputs: true}).combine({
              reducer2: ee.Reducer.stdDev(), sharedInputs: true}).combine({
              reducer2: ee.Reducer.mode(), sharedInputs: true}).combine({
              reducer2: ee.Reducer.median(), sharedInputs: true}).combine({
              reducer2: ee.Reducer.mean(), sharedInputs: true});

print(ui.Chart.image.series(collection, KansasMap, reducers, 200));


texts: Applying https to Carto map but requesting insecure data?

texts:
I am building a site using data from CartoDB that locates someone on a map, and Chrome requires that be over https. I added an SSL, but it was not working, the page showed up black. I found out that it was calling things that were not over https, and fixed a couple of them by adding them to my site, but I still get the error below. Now it locates someone, but does not show the map data.

The data should show up like the line here (permission denied on location, because not https):
http://bakkenpipelinemap.com/

But https shows the location, but not the data:
https://bakkenpipelinemap.com/

The error below from the Console says it's requesting something, looks like it's the map data, which is over http not https. Any ideas?

  Mixed Content: The page at 'https://bakkenpipelinemap.com/' was loaded
  over HTTPS, but requested an insecure script
  'http://nittyjee.carto.com/api/v1/map?stat_tag=6c459170-8476-11e6-af6f-0e05a…tivity%22%3A%5B%22cartodb_id%22%5D%7D%7D%5D%7D&callback=_cdbc_2698715151_1'.
  This request has been blocked; the content must be served over HTTPS.


texts: How to register own CQL filter in GeoServer?

texts:
I have created new FilterFactory implementation, added it to META-INF/services in file org.opengis.filter.FilterFactory.
Added Interface extending Filter.
Added implementation of this Interface.
When I try to call it in cql_filter I get CQLException "Function not found"..
What am I doing wrong?


texts: Buffering around raster using gdal and numpy?

texts:
I have a stream raster with stream cells containing stream-order values. I need a buffer of 1000m around the streams and it would be good if the values in the buffer raster contain the stream-order to which it is close to.

I am using gdal and numpy modules to handle raster.

How can I do this?


texts:
I was hoping someone could help me with this. I know this is probably easy but I am still new to ArcMap.

I need to label pipe mains with the year installed but my attributes have the entire date.  How can I format the label to only display the year (yyyy)? I am using ArcMap 10.2.2 and VBScript. I have included screenshots to show where I am at currently.


texts: Get the value/description from a domain or subtype for a given field's value?

texts:
I need to get a value from an IObject.Fields.Field using ArcObjects 10.x using the criteria listed below:

If the field has a subtype associated with it, I need to get the
subtype name/description based on the field value.
If the field has a coded value domain associated with it, I need to
get the domain's name/description based on the field value.
If the field has neither a subtype or domain associated with it, I
want the field value itself returned.
Having this functionality wrapped up in one call like "ReadFieldValue(IObject obj, int fieldIndex) is preferred.


texts: Intersecting multipolygon layer with polygon in GeoDjango

texts:
I have a Region model (deriving from django.contrib.gis.db.models.Model), with a geom field, which is a django.contrib.gis.db.models.MultiPolygonField.

I'd like to compute the intersection of this model (representing a set of polygonal regions stored in PostGIS) with a GeoJSON polygonal query (in the context of a web service):

from django.contrib.gis.geos import GEOSGeometry

query = GEOSGeometry(
    '{"type":"Polygon","coordinates":[[[-80.983786,43.929011],[-80.511513,43.778458],[-80.291852,44.079184],[-80.775108,44.232127],[-80.983786,43.929011]]]}'
)

I tried many things:

results = Region.objects.all().intersection(query)

or

from django.contrib.gis.db.models.functions import Intersection

results = Intersection(Region.objects.all(), query)

but it seems I don't have the right approach.

Note that I do not want to simply compute the subset of regions which intersect with the query, but rather their exact intersection (in other words the resulting set of polygons should have the exact same bounding box as the query).


texts: How do I set the expiration of tiles globally in Geoserver? In the tilecache tab for a layer it says setting the value to 0 (default) uses the source setting. Where is this source setting set? I have tried altering geowebcache.xml but this is just per layer basis.


texts: Performing join on simultaneous spatial and attribute criteria in ArcGIS for Desktop?

texts:
I have two point layers and want to perform a spatial join on nearest distance.

Normally it is easy with ArcMap. But now I have a restricting condition on the spatial join. The restriction only allows the attributes to be joined if the nearest point is within the same country.

Both point layers have information about the country in their attribute table (e.g., "de" for Germany). So it must be some spatial join like:

Join attribute B to layer A on nearest distance where countryID_Layer A = countryID_Layer B.

Does anybody have an idea how I could solve this?

EDIT: Image I have following attributes: Table1(hotels): h_key, l_iso_a2, geom and in Table2(iata_codes): loc, ctry, geom

EDIT2: I tried to write the SQL-Statement, but it does not work somehow. Can anybody help me with this statement?:

select distinct on (h.h_key) h.h_key, h.l_iso_a2, i.loc, i.ctry, distance
from ( select h.h_key, h.l_iso_a2, i.loc, i.ctry as country2, st_distance (h.geom,
i.geom) as distance
from hotels AS h, iata_codes AS i
where h.l_iso_a2 = i.ctry order by h.geom <-> i.geom ) as iata_codes_h_key;


texts: Sum area of multiple polygons of same layer that are all contained within different polygon

texts:
I have multiple polygons of the same layer that are all contained within another polygon of a different layer (please see the example below- picture is of Great Barrier Island and surrounding islands, in New Zealand).

I am hoping to find a way to sum the area of the pink polygons (government owned land) that are contained within the larger green polygons (the actual island projection) in ArcMap. As can be seen in the figure, there are multiple islands, each of which contain varying amounts of government owned land.

Ideally, I am looking for a way to calculate this information for many islands within the New Zealand archipelago. If there is a way to calculate this en masse, that'd be great, but if it has to be done manually, that's okay too.

I am relatively new to GIS and ArcMap. Also, I am proficient in R (but not Python) if that can assist the analysing process.


texts: Access Stamen tiles flavors with OpenLayers

texts:
I'm using OpenLayers (v2.13.1) and would like to use Stamen Toner tiles as background.

The example is simple to understand: http://maps.stamen.com/test/openlayers.html

But I'd like to use the "lite" variant of the toner tiles, how is it possible with OpenLayers?