texts: SQL query within Select tool in ArcGIS model going very slow

texts:
I have been working on finding the shortest routes from a series of residential blocks (a shapefile with polygon centroids) to a a number of recreational areas. I have a total of 9 resulting OD cost matrices (Playing fields within 1000m, Playing fields within 600m, Playing fields within 300m, and then  recreational areas and green recreational areas all over again).

Since each destination polygon (playing fields, gren rec areas, rec areas) has multiple points of entry, the OD cost matrix contains multiple routes to each destination. I only need the shortest of those routes.

I have therefore built a model to go through each of the OD cost matrices and select the shortest route from each block to each listed destination.

I've only done about three models before this, and I feel it is quite crude.

The model is comprised of the following:

Part 1 - Iterates through the OD Cost Matrix files
Part 2 - Nested in Part 1 -Exports all routes from a single block to a separate shapefile using Iterate row selection and Copy feature
Part 3 - nested in Part 2 - Here the model iterated through the destination ID's (each destination polygon has its own id, as do all its perimeter destination points) , selects all entries with the same destination ID's and saves them to a new shapefile. Then it takes that shapefile, uses 'Calculate value' to select the smallest distance value, and using another Select tool with the following SQL query selects the entry with the smallest distance and saves it as a separate shapefile.

"Total_Leng" = %output_value%

where %output_value% is the value passed by the 'Calculate value' tool

This last selection tool takes from 3 to 6 minutes per iteration. The model has been running for 6 days now, and has only gone through 600 of 2400 entries from the FIRST of nine shapefiles.

The biggest problem is Part 3
I am sure there is a more efficient alternative to saving selections as shapefiles. I would think there would be an option to keep selections in memory, but I really cannot find it and I read that model builder isn't so great at working from memory.

Any suggestions as to how I can improve performance? The results of this model are good, but it is simply too slow. I have been searching for similar problems but I haven't found anything specific. Maybe someone here has any suggestions? I know I could do a better job with Python but I cannot afford the time to get into Python beyond the basics (which is what I know - the very basics).

I've added images of the model below, I could also add the python scripts if required.


texts: How to "clean up" undesired values in a slope raster

texts:
I have a reclassified slope raster with 3 values in band 1. My problem is that I want to polygonize this layer, but the surface is so irregular with a lot of plots in terraces and because of this it becomes in a vector layer with infinite polygons.

What can I do to clean up the raster layer to obtain bigger polygons of slope?


texts: How to create "rotate canvas" shortcuts in QGIS?

texts:
Is it possible to create "rotate canvas" shortcuts in QGIS?
I'm thinking about something like "Ctrl + mouse scroll up" to rotate the canvas clockwise and "Ctrl + mouse scroll down" to rotate the canvas counter clockwise...


texts: Trying to automatically update, getting empty output instead

texts:
I am trying to make a model that basically updates an original file with new data. I do this by using a mask to "erase" old data and then "merging" in the new data, however this creates a separate feature class. Since the original file is on a shared network, I can't just overwrite the file.

To achieve my goal of updating the original file, I empty it with "Truncate table"(Input %Name% being the orginal file, decided in the "Iterate Feature Classes"), after which I use "Append" to put the new data inside.

Seeing as my users might want to run multiple feature classes at once, I use an "Iterate Feature Classes" at the beginning of my tool.

My problem is in the result, my "updated feature class" gets created correctly with data and, for as far as I can tell, the "original feature class" gets emptied. However, somehow when the tool is done running, it returns an empty output (over the entire feature class, not just the updated bit underneath the mask). I have tried all different options of preconditioning, but nothing seems to be a solution.

It is important to note that, if I run these bit by bit in ModelBuilder itself, it does produce the correct result, but when I run the model it does not.

The same result happens if I use "delete features" or "delete rows"
Putting them into seperate models (making the updated file, emptying the original file and then appending) and putting those models into one model, also gives the same empty output

I am running ArcGIS Desktop 10.3.1.
As extra information, I added a pictures of the attribute tables before and after the model runs and the model itself. I have been at this for several days and I cannot figure it out. The inout "New file" is the adjustments made feature class from the screenshot below.

Before

Adjustments made

After

Model


texts: Converting Map Coordinates to Layout Coordinates in rotated Dataframe?

texts:
I am trying to convert Map X and Y coordinates to corresponding Layout coordinates in a rotated dataframe.

The code found in this answer to Convert point XY to page units XY using arcpy? works fine, but only for unrotated dataframes.

I've also found a modification of the above code in this post on Reddit, but that doesn't seem to work either. I believe the code was making use of a rotation matrix as discussed in this post, but without success, in my opinion.

The results I get are somewhat close when the dataframe is not rotated very much, but they exaggerate the greater the rotation.

Has anyone had definitive success with converting X and Y coordinates to Layout coordinates?

I'm looking for a python solution within ArcGIS. I'm using 10.3.1 Standard License.


texts:
I am trying to automate a map I create every so often. The map flags counties based on certain commodities, and color codes them based on price changes. I am able to get each layer created, but would like to reorder the commodities so the "red" commodities are on top, orange below that, low risk below that, and then move the State layer to the top.

My code is below - the error result I get is:

Assertion Error: Did not find move layer

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

arcpy.AddJoin_management("County_Projected_AEAC", "ID", "C:/Users/SS/Desktop/ArcGIS Maps/MTW/Spring_2020.xlsx/Spring_2020_MTW$", "modified_FIPS", "KEEP_ALL")
arcpy.Select_analysis("County_Projected_AEAC", "C:/Users/SS/Desktop/ArcGIS Maps/MTW/Beef_Cattle_Python.shp", '"Spring_2020_MTW$.Flag_Beef_Cattle" = 1')
arcpy.Select_analysis("County_Projected_AEAC", "C:/Users/SS/Desktop/ArcGIS Maps/MTW/Feed_Cattle_Python.shp", '"Spring_2020_MTW$.Flag_Feed_Cattle" = 1')
arcpy.Select_analysis("County_Projected_AEAC", "C:/Users/SS/Desktop/ArcGIS Maps/MTW/Corn_Python.shp", '"Spring_2020_MTW$.Flag_Corn" = 1')
arcpy.Select_analysis("County_Projected_AEAC", "C:/Users/SS/Desktop/ArcGIS Maps/MTW/Cotton_Python.shp", '"Spring_2020_MTW$.Flag_Cotton" = 1')
arcpy.Select_analysis("County_Projected_AEAC", "C:/Users/SS/Desktop/ArcGIS Maps/MTW/Dairy_Python.shp", '"Spring_2020_MTW$.Flag_Dairy" = 1')
arcpy.Select_analysis("County_Projected_AEAC", "C:/Users/SS/Desktop/ArcGIS Maps/MTW/Hogs_Python.shp", '"Spring_2020_MTW$.Flag_Hogs" = 1')
arcpy.Select_analysis("County_Projected_AEAC", "C:/Users/SS/Desktop/ArcGIS Maps/MTW/Poultry_Python.shp", '"Spring_2020_MTW$.Flag_Poultry" = 1')
arcpy.Select_analysis("County_Projected_AEAC", "C:/Users/SS/Desktop/ArcGIS Maps/MTW/Soybeans_Python.shp", '"Spring_2020_MTW$.Flag_Soybeans" = 1')
arcpy.Select_analysis("County_Projected_AEAC", "C:/Users/SS/Desktop/ArcGIS Maps/MTW/Wheat_Python.shp", '"Spring_2020_MTW$.Flag_Wheat" = 1')

States = arcpy.mapping.ListLayers(mxd)[-1]
County = arcpy.mapping.ListLayers(mxd)[-2]
Beef_Cattle = arcpy.mapping.ListLayers(mxd)[-3]
Feed_Cattle = arcpy.mapping.ListLayers(mxd)[-4]
Corn = arcpy.mapping.ListLayers(mxd)[-5]
Cotton = arcpy.mapping.ListLayers(mxd)[-6]
Dairy = arcpy.mapping.ListLayers(mxd)[-7]
Hogs = arcpy.mapping.ListLayers(mxd)[-8]
Poultry = arcpy.mapping.ListLayers(mxd)[-9]
Soybeans = arcpy.mapping.ListLayers(mxd)[-10]
Wheat = arcpy.mapping.ListLayers(mxd)[-11]

low_risk_commodities = [Poultry]
for x in low_risk_commodities:
    x.visible = False
    arcpy.RefreshTOC()

orange_commodities = [Wheat, Hogs, Dairy, Cotton, Corn, Feed_Cattle, Beef_Cattle]

for x in orange_commodities:
    arcpy.mapping.MoveLayer(df, arcpy.mapping.ListLayers(mxd)[0], x)

The code is correctly turning off the poultry layer but will not move the layers due to the error. I know it is likely in how I'm defining the layers, but I'm not sure how to properly do this.

When I type in Wheat, it gives me a little purplish icon, indicating it is maybe not a layer but is maybe a shapefile?


texts:
<Rule>
      <Title>EU-SPA</Title>
      <ogc:Filter>
      <PropertyIsEqualTo>
         <PropertyName>BEWERT2006</PropertyName>
         <Literal>EU-SPA</Literal>
      </PropertyIsEqualTo>
      </ogc:Filter>
      <PolygonSymbolizer>
        <Fill>
          <GraphicFill><Graphic>
            <Mark>
              <WellKnownName>shape://slash</WellKnownName>
              <Fill>
                <CssParameter name="fill">#66ff33</CssParameter>
                <CssParameter name="fill-opacity">.7</CssParameter>
              </Fill>
              <Stroke>
                <CssParameter name="stroke">#000000</CssParameter>
                <CssParameter name="stroke-width">1</CssParameter>
              </Stroke>
              <Stroke />
            </Mark><Size>5</Size>
            </Graphic></GraphicFill>
        </Fill>
        <Stroke/>
      </PolygonSymbolizer>
    </Rule>

I would like to have a square with slanted lines inside in which the fill must be one colour and the stroke must be another. I have written this code to make it run but it just shows without any colour. Can somebody share some ideas?


texts: Leaflet js style geojson layer according to zoom level

texts:
Leaflet js beginner following/adapting the choropleth example. I need to get something like this mapbox map where the feature styles change according to zoom level. I am loading a line geojson that gets style from getColor / getWidth functions related to feature.properties. How do I get the weight to be 1 when currentZoom is <= 16?

var geojson = L.geoJson(roads, {
      style: style
    }).addTo(map);

function getColor(p) {
  return p >= 100 ? '#00F7BD' :
    p >= 40 ? '#24abc9' :
    p >= 20 ? '#ef561a' :
    p >= 0 ? '#4D2A2D' :
    '#feffbd';
}

function getWidth(w) {
   return w > 20 ? 4.8 :
          w > 17 ? 4.5 :
          w > 16 ? 4 :
          w > 15 ? 3.5 :
          w > 14 ? 2.5 :
          w > 10 ? 1.5 :
          w > 0 ? 0.5 :
          0;
}

// map.on('zoomend', function () {
  var currentZoom = map.getZoom();

function style(feature) {
      if (currentZoom <= 16) {
        return {
          color: getColor(feature.properties.paP),
          weight: 1,
          opacity: 0.5,
        };
      }
      else {
        return {
          weight: getWidth(feature.properties.stW),
          color: getColor(feature.properties.paP),
        };
      }
}
// });

Second approach inspired by Stephen's answer @Stephen Lead

   function style(feature) {
      var fa;
      return {
          color: getColor(feature.properties.paP),
          weight: fa * getWidth(feature.properties.stW),
          opacity: 0.6,
          lineCap: 'square',
          lineJoin: 'square',
        };
    }

    map.on('zoomend', function () {
        var currentZoom = map.getZoom();
            if (currentZoom <= 16) {
                  fa = 0.5;
            }
            else {
                  fa = 3;
            }
    });


texts: Autodesk MapGuide Enterprise MrSid on Linux

texts:
I am running MapGuide Enterprise 2011 on CentOS 5 Linux. Using the Autodesk Raster FDO, I loaded a mrsid file (through the "Load file-based data" method). Whenever i try to view the layer with the mrsid file, the server crashes. most of the time it says "Segmentation Fault", but sometimes it just dies without saying anything.

Here's what the error log says when I try to view the layer:

<2011-01-14T10:39:25>   53263248        Autodesk MapGuide Studio v2.2.0.5305       198.245.206.237 Administrator
Error: An exception occurred in FDO component.
    Only MOSAIC(), CLIP() and RESAMPLE() are supported.
StackTrace:
- MgFeatureServiceHandler.ProcessOperation() line 83 file FeatureServiceHandler.cpp
- MgOpSelectFeaturesSpatial.Execute() line 110 file OpSelectFeaturesSpatial.cpp
- MgServerSelectFeatures.SelectFeatures() line 266 file ServerSelectFeatures.cpp

I've tried a few other mrsid files with the same problem. TIFF files load fine.

Anyone has any insight on this?


texts: I have a problem when trying to load raster data (a TIF) in ArcMap. ArcMap 10 instantly shuts down and the well known window appears "ArcGis Desktop has encountered a serious application error and is unable to continue.Does anyone know why? What is the solution? i have windows 7 32 bit,1giga video card .


texts: Is there a way to control label density in mapserver?

texts:
I'm reading a US County Tiger Shapefile and everything is plotting correctly, but I want to slim down the number of labels shown at scales larger than the state level because it begins to overwhelm the map.  Is there a way to modify this to do that, but leave the zoomed in look and feel alone?

LAYER
NAME 'us_counties'
TYPE POLYGON
DUMP true
TEMPLATE fooOnlyForWMSGetFeatureInfo
DATA './shapefiles/uscounties_20130327.shp'
METADATA
  'ows_title' 'counties_us'
END
PROCESSING "LABEL_NO_CLIP=True"

PROJECTION
  'proj=longlat'
  'ellps=GRS80'
  'datum=NAD83'
  'no_defs'
END

LABELITEM 'NAME'
CLASS
  NAME 'uscounties_20130327'
  STYLE
    WIDTH 0.91
    OUTLINECOLOR 0 0 0
    #COLOR 108 180 98
  END
  LABEL
  FONT verabd
  TYPE truetype
  SIZE 12
  COLOR 0 0 0
  OUTLINECOLOR 255 255 255
  ANGLE 0
  POSITION cc
  FORCE true
  ANTIALIAS true
  PARTIALS false
END
END
  END


texts: Automatically generating label lines in QGIS?

texts:
I've been experimenting with orthographic projections and I've run into a problem with the labels.

I would like to label the countries but I want the labels to run west-east along the country which isn't (currently) an option in the QGis labeling options.

So, I know I can place them along a line but I lack the patience to draw a line in each country. So I'm looking for an automatic way to draw a line through the centroid of a polygon and to the edges so that I can use it as a label placement line.

By 'label placement line' I mean I need a line that crosses the country that I can use as the basis of my labels so that I can use follow line (parallel I think) in the labeling panel. (see for example this question). So they will be a straight horizontal line in EPSG:4326 but nicely curved in my orthographic projection (+proj=ortho +lat_0=51.470129 +lon_0=-0.452751 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs).

Ideally it would be a dense line so it curves when reprojected but I can always run a densify process on the output.


texts: Creating flowline by connecting points in ArcGIS Desktop?

texts:
I have a dataset that contains X,Y coordinates of a glacier and the surface elevation is plotted as a Z-coordinate (see attachment).

What I want to do is connect the top with the bottom of the glacier by a line, that goes perpendicular to the contours (a flowline). Does anyone know how this can be done? I thus want to create a transect (line) perpendicular to the contourlines...


texts: I want to use ArcGIS for viewshed analysis on a Mac. Would virtual machine software such as Parallels or VMware Fusion be powerful enough for viewshed analysis compared with bootcamp?


texts:
Is there any means to combine the function of zoom to point and place a pin plugins in QGIS? Would it be possible to enter a xy co-ordinate and a pinpoint is placed at the point of the grid reference. My apologies if I have missed something obvious on other forums. GSherman maintains these both so I am hoping this query reaches him!?

Thanks


texts: Downloading and merging 0.5 minute WorldClim data for (long) list of coordinates using R

texts:
I want to download temperature and precipitation data from WorldClim for a huge list of coordinates.  I am working in R, and have a solution using getData() from raster(), but it's going to take forever for the giant list of coordinates I have.  I posted a similar question over on stackoverflow in hopes I could streamline the looping.  However, I think even if I can streamline the looping, this might take forever.
Does anyone have a suggestion for a more efficient way of doing this?
I've Googled a whole bunch, and can't find any solutions in R.  I do have a reproducible working example:
# Download and merge 0.5 minute MAT/MAP data from WorldClim for a list of lon/lat coordinates
# This is based on https://emilypiche.github.io/BIO381/raster.html

# Make a dataframe with coordinates
coords <- data.frame(Lon = c(-83.63, 149.12), Lat=c(10.39,-35.31))

# Load package
library(raster)

# Make an empty dataframe for dumping data into
coords3 <- data.frame(Lon=integer(), Lat=integer(), MAT_10=integer(), MAP_MM=integer())

# Get WorldClim data for all the coordinates, and dump into coords 3
for(i in seq_along(coords$Lon)) {
  r <- getData("worldclim", var="bio", res=0.5, lon=coords[i,1], lat=coords[i,2]) # Download the tile containing the lat/lon
  r <- r[[c(1,12)]] # Reduce the layers in the RasterStack to just the variables we want to look at (MAT*10 and MAP_mm)
  names(r) <- c("MAT_10", "MAP_mm") # Rename the columns to something intelligible
  points <- SpatialPoints(na.omit(coords[i,1:2]), proj4string = r@crs) #give lon,lat to SpatialPoints
  values <- extract(r,points)
  coords2 <- cbind.data.frame(coords[i,1:2],values)
  coords3 <- rbind(coords3, coords2)
}

# Convert MAT*10 from WorldClim into MAT in Celcius
coords3$MAT_C <- coords3$MAT_10/10


texts: Dissolving Multiple Ring Buffer by Buffer Ring Iteration Order using ArcPy?

texts:
I have a lines feature class with the following attribute table:

When I run the MultipleRingBuffer_analysis function using the following code (provided by user @Emil Brundage), buffers are drawn using the sum of the "width" column of the attribute table and a list of distances 10, 20 and 30 meters.

import arcpy
arcpy.env.workspace = "C:/data/Py.gdb"
fc = "lines"
distances = [10,20,30]
sr = arcpy.Describe(fc).spatialReference
#create feature class
arcpy.CreateFeatureclass_management("C:/data/Py.gdb",
                                     "Multibuffers",
                                     "POLYGON",
                                     spatial_reference = sr)
#insert cursor to insert geometries into feature class
with arcpy.da.InsertCursor ("Multibuffers", ["SHAPE@"]) as iCurs:
    with arcpy.da.SearchCursor(fc, ['SHAPE@', 'width','OID@']) as cursor:
        for geom, width, oid in cursor:
            #buffer in memory
            arcpy.MultipleRingBuffer_analysis(Input_Features=geom,
                                             Output_Feature_class=r'in_memory\buffs',
                                             Distances=[i+width for i in distances],
                                             Buffer_Unit='meters')

            #iterate buffers and add to output feature class
            with arcpy.da.SearchCursor (r'in_memory\buffs', ["SHAPE@"]) as sCurs:
                for row in sCurs:
                    iCurs.insertRow (row)

            #delete in memory buffers
            arcpy.Delete_management (r'in_memory\buffs')

Output is one non dissolved feature class containing all the buffers as POLYGON objects (screenshot 1). I need to find a way to dissolve the buffers by ring iteration order in order to obtain something similar as screenshot 2 but with the right distances applied.


texts: Orthophoto processed with Photoscan distorted in QGis

texts:
I processed an orthophotomosaic via Photoscan using geotagged images in EPSG 4326 and it is displayed correctly in Photoscan and GoogleEarth. When I put it in QGis it is distorted. The coordinate system in QGis is set correctly.

Does anyone know what could be the issue?


texts: How do I select objects I click on in ArcGIS (.Net)

texts:
I have a MapView in my xaml code, which contains several GraphicsLayer. The first layer is the actual map, and the second is small SimpleMarkerSymbol representing locations.

What do I do to be able to click on them?

I use the trigger MapViewTapped which gives me MapViewInputEventArgs with coordinates of my click, and I have my MapView. According to the developer guide, I should be able to do the following:

var objectIDs = await lyr.HitTestAsync(MyMapView, myClick.Position);

lyr.SelectGraphics;

However, nothing is found. Note: the link uses FeatureLayer instead of GraphicsLayer.
Unlike FeatureLayer, GraphicsLayer does not have a setter in its SelectGraphics property.

So - how do I click on an marker and have it returned?


texts: Connection between ArcGIS Desktop and ArcGIS Online

texts:
When connecting ArcMap with ArcGIS Online by importing the map from ArcGIS Online in ArcMap on my desktop and using 'Synchronize Local Edits with Server' everything that I edit on ArcMap is also edited online. But, when I try to add another column to my attribute table I cannot see this change in ArcGIS Online (although I do refresh and save etc). Also the other way around, when I add another column in the data of the map in ArcGIS online, I cannot see this in ArcMap.

Does anybody have a solution for this?


texts: How to duplicate a map layer at most once as the bottom layer in qgis3 (python)

texts:
I wanted to duplicate a layer in qgis3 using python (display the same vector data as different layers) per this question (How can I duplicate a layer in QGIS?). In addition, I need to:

make sure that the added clone layer is at the bottom of the layers
check if the clone already exists and duplicate only if the clone doesn't exist among the current list of layers (to avoid multiple clones).

The code I have so far (for duplicating the first layer) is:

vl = iface.mapCanvas().layer(0)
vl_clone = iface.addVectorLayer(vl.source(), vl.name() + "_clone", vl.providerType())

Currently the new layer is added to the top and becomes layer(0). I tried to follow existing GIS.SE questions (e.g. How to move layers in the Layer Order Panel using PyQGIS?). But existing questions/answers are for QGIS 2, and involve breaking changes in QGIS3 (e.g. iface.layerTreeCanvasBridge())

Does anyone know how to extend the above code to work in QGIS3?

-- EDIT --

The current code I have so far from synthesizing existing answers/documentation is below. Still, the layer order doesn't change. I tried different values in insert() such as 3, to no avail.

lyr1 = iface.mapCanvas().layer(0)
lyr1_clone = iface.addVectorLayer(lyr1.source(), lyr1.name() + "_clone", lyr1.providerType())

bridge = iface.layerTreeCanvasBridge()
treeRoot = QgsProject.instance().layerTreeRoot()
treeRoot.setHasCustomLayerOrder(True)
order = treeRoot.customLayerOrder()
order.insert( -1 , order.pop() ) # lyr1_clone to the bottom # lyr1_clone to the bottom
treeRoot.setCustomLayerOrder( order )
bridge.setCanvasLayers()


texts: The maps look good, and the legend symbols match the maps.  Once converted to .pdf the map legend symbols degrade to basically nothing.  Horizontal red line patterns, or any colored line patterns for that matter are degraded to "pink" or "colored" legend boxes after conversion to PDF, which makes the legend useless.


texts: Two vector layers, same CRS, wrong coordinates

texts:
I have two vector shapes, downloaded from government providers, that should be overlapping perfectly but one appears way of to the east of the other one. They are both set to the same CRS.

I work with several other layers too, which all fit perfectly with the first one.  When I load in the second layer (the one that is off) it says "CRS was undefined". I then set it to the project CRS manually.

Any idea what I'm doing wrong here?


texts: I want to draw a closed line in QGIS1.8,but I don't know how to make the last point coincide with the first point.What am I supposed to do?


texts: How to add select/deselect all layers button on the leaflet webmap created via qgis2web plugin?

texts: I have a leaflet web map created via qgis2web plugin using qgis software. There is no option in plugin window to generate select/deselect button for the map layer. Now I want to add that button on the map, but I don't know which code part and where I should add it.


texts: How to calculate Cost Surface Analysis in ArcGIS 10.x

texts: I have about 10 points and I want to run Cost Surface Analysis. I’m trying to find what distance a person can walk in a given amount of time from a certain point (5-10 or 15 minutes). What procedure follow and what data I need. I have DΕM of the region. I use ArcGIS 10.3 and I read from online help the Creating a cost surface raster, but I'm not sure what I have to do.


texts: Why do Print composer templates not save my maps in their entirety?

texts: I have been using QGIS and its various updates in a limited way to make maps for a few years now, but I have always had the same problem with the Print Composer.  I make my layered maps with grouped annotations and save them as templates, which I'd always assumed would allow me to recall them at another time, when the map would pop up just as I'd left it and if I wanted to make minor changes or adjustments, I could then do so. This has never happened. What does happen is that when I open Print Composer files, the annotations appear, but not the original map images. I have looked at various QGIS manuals but cannot find the answer to this. Would someone be able to explain to me if there is a way to save the complete content of maps for retrieval and/or editing at a future date?  To me, this would appear to be the main function of a GIS, so I must be misunderstanding how it works.


texts:
I have flown topo lines covering a large AOI. Some of these lines stretch across the whole area and therefore aren't closed while others just cut off and also aren't closed. I want to calculate the area of the lowest 5 contour lines per each basin (polygon).

I tried extending the lines and then converting them into polygons but this didn't work because some of them did not close.

How do I find the lowest 5 lines per each basin (there are about 35)?

I am using ArcMap 10.3 and AutoCAD Civil 3D. My company does not have the spatial analyst package.


texts: Find feature classes that are not in one database and copy them to another

texts:
I have a database that meets a particular data standard and another that only partially meets this same standard. I am trying to combine data from the partially standardized database to the standardized database without creating repeats.

Using the Python sets module seems like it may provide a solution. The code below only provides me with an empty list or a list that doesn't contain the intended feature classes. Once I get this code to work, my intention is to use the produced list to insert into something similar to arcpy.FeatureClassToFeatureClass_conversion (migrate, r'C:\database2.gdb', filename)

import arcpy, os
from arcpy import env

# List data that is in DB1 (partially compliant DB)
workspace = r'C:\database1.gdb'
FCs_in_DB1 = []
for dirpath, dirnames, filenames in arcpy.da.Walk(workspace,
                                                  datatype="FeatureClass",
                                                  type="Any"):
    for filename in filenames:
        FCs_in_DB1.append(os.path.join(dirpath, filename))

# List data that is in DB2 (compliant DB)
workspace = r'C:\database2.gdb'
FCs_in_DB2 = []
for dirpath, dirnames, filenames in arcpy.da.Walk(workspace,
                                                  datatype="FeatureClass",
                                                  type="Any"):
    for filename in filenames:
        FCs_in_DB2.append(os.path.join(dirpath, filename))

# What I have tried so far.
# I have also switched the 1 and 2 in list names in attempt to get desired results.
#migrate = list(set(FCs_in_DB2) - set(FCs_in_DB1))
#migrate = [x for x in FCs_in_DB2 if x not in FCs_in_DB1]
#migrate = FCs_in_DB1.symmetric_difference(FCs_in_DB2))
#migrate = list(set(FCs_in_DB1) ^ set(FCs_in_DB2))

print migrate

I am hoping to find what is wrong with my code or perhaps be led to a better approach to this problem altogether?


texts: Scale difference between WGS84 WebMercator and Lambert projection

texts:
There's something I don't understand: I have a map in ArcMap in Lambert RGF93 projection (conic projection system for France) at scale 1:10000 then I change the coordinate system to WGS84 WebMercator map, I get a map at 1:10000 but data are zoomed in.

Here's my map in RGF93/Lambert-93:

And here's the WebMercator one:

I understand it can be differences in angles or areas between both projections but why the same scale doesn't give similar zoom level?


texts: Creating layers depending on a property?

texts:
I have a .js file with coordinates for internships:

var internships = [{
  "features": [
    {"type":"Feature","properties":{"category":"entretient","Name":"green"},"geometry":{"type":"Point","coordinates":[50.807149, 3.162994]}},
    {"type":"Feature","properties":{"category":"securité","Name":"blue"},"geometry":{"type":"Point","coordinates":[50.334421, 3.290146]}},
    {"type":"Feature","properties":{"category":"secretaria","Name":"red"},"geometry":{"type":"Point","coordinates":[50.744787, 2.256216]}}
  ]
}];

I've found this bit of code allowing me to create layers depending on a property and here what my JS looks like:

$.getScript("CoordinatesPdC.js");

function mapLoad() {
  var sécuritéLayer = new L.LayerGroup();
  var secrétariatLayer = new L.LayerGroup();
  var entretientLayer = new L.LayerGroup();

  var map = L.map('map').setView([50.2910, 2.7775], 8);

  L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    maxZoom: 18,
    attribution: 'Map data &copy; OpenStreetMap contributors, ' +
    'CC-BY-SA, '
  }).addTo(map);

  var marker = L.marker([50.2910, 2.7775]).addTo(map);

  var entretientLayer = L.geoJson(internships, {
    filter: function (feature, layer) {
      return (feature.properties.category === "entretient");
    }
  }).addTo(map);

  var sécuritéLayer = L.geoJson(internships, {
    filter: function (feature, layer) {
      return (feature.properties.category === "sécurité");
    }
  }).addTo(map);

  var secrétariatLayer = L.geoJson(internships, {
    filter: function (feature, layer) {
      return (feature.properties.category === "secrétariat");
    }
  }).addTo(map);

}

window.onload = mapLoad;

But now I have to create the markes assigned to these layers, how can I achieve that?


texts: How to resize an infoWindow dynamically, based on the contents

texts:
How can I set the size of an infoWindow dynamically, based on the length of its Title, when the Title isn't known until the user clicks on a map feature?

Take the map here as an example - how could you resize the popup infoWindow to accommodate varying Title lengths?

This map uses the ${Name} notation to represent the Title field, with the field value substituted into the infoWindow when the feature is clicked.

The infoWindow has a fixed size (set using map.infoWindow.resize(250,95)). Since the length of the title isn't known until the feature is clicked, how is it possible to resize the infoWindow based on the length of the Title?

Using ArcGIS Server 10 with the JavaScript API 2.3


texts: Reprojecting in postGIS from 900913 to 4326. Wrong coordinates

texts:
I'm working with postGIS 2.1
I have a table called 'xacementos' created with the srid: 900913. I took the points coordinates from google maps and I want to convert them to SRID:4326 to show the points in a map created over leaflet. (if a input the 900913 coordinates as geoJSON in the leaflet code the points show a little bit displaced)
I'm using the function ST_Transform:

SELECT ST_AsText(ST_Transform(geom,4326))
FROM xacementos
WHERE enderezo='Rua Franja 9-11';

but the coordinates that show up are not correct:

"POINT(-7.54315344075162e-005 0.000389606525237254)"

The input coordinates in the 900913 srid are "POINT(-8.397 43.3708)"

I'm newby in this postGIS stuff and I'm not sure if i'm following the right path.

Any suggestion??


texts: Anyone have any idea how I could connect the GNIS list of waterbodies to StationID In the NOAA dataset?

texts:
We are building lookups of waterways across the country. Have built a lot of info and trying to automate links to stationid in the gage system of NOAA.  Does anyone know where to get the waterbody id? It seems to be a text field.
We cannot find a dataset anywhere to link the two.


texts: Non-Earth maps, tile server and scale in leafletjs?

texts:
I am looking at putting together a tile server for other astral bodies, such as Mars and the moon. At the moment I am using leafletjs and the Control.Scale at the client end (though could be any client API), which assumes Earth as the planet.

Should the planet radius be provided by the server to the client or should this be specified in the client? If it is the server, is there a metadata file that would be looked for and what is its format? If it is the client, should this be done via a custom TileLayer?

At the moment I am looking to use an open street maps style URL, such that:

http://{s}.tile.myhost.mydomain/planet/<planet_name>/<layer>/{z}/{x}/{y}.png


texts: From many raster tiles to one geopackage

texts:
I need to unify many raster tiles of global land use (856 tiles, in a total of 297 GB, single band, 8 bite unsigned integer, WGS84) in a single file, I chose geopackage because of the high dimensions.

I tried using the following steps:

building a GDAL virtual raster with all the tiles (Parameters: NOT checking the separate band option, high resolution, resempling=mode)
converting it into a geopackage by means of gdal_translate:

gdal_translate -of GPKG -co QUALITY=100 -b 1 -mask 1 -r mode -a_offset value input.vrt output.gpkg

but the output has the following errors:

there are some pixels with a wrong value
the output raster has 4 bands (3 with the same values and the last with NoData, 255)

I noted that in every passage there is a shift of the pixels, a sort of cell value resampling, despite mode resampling should select the value which appears most often of all the sampled points.

How it is possible to create a single file with exactly the same pixels and parameter as the input rasters?

The original raster tile:

the same tile converted to virtual raster (info on the same cell):

the final GPKG output (info on the same cell):

In the GPKG there is a wrong new cell value (for ex. 41 instead of 40) and 4 bands instead of 1.

Note that despite the partial wrong output (some wrong pixel values and presence of 4 bands instead of 1), from a dataset of 297 GB I obtained a geopackage of 10 GB for the southern hemisphere and 34 GB for the northern hemisphere!


texts: Cursor not centered on feature - OpenLayers

texts: Can anyone tell me why the cursor is not centered on the feature? This can be seen in the Overlay event of the map and even in the onClick event. We can only click on the element in a certain position, never in the center.


texts: Loading data from shapefiles into specific feature classes

texts: I am trying to bulk load a collection of shapefiles into empty feature classes in ArcMap 10.2.  The issue here is that the names of the shapefiles and the feature classes need to match up.  Is there any way to make this happen in Model Builder or arcpy?


texts: I am relatively inexperienced in using GIS systems, but was able to make great use of QGIS several years ago to make distribution maps for my dissertation. I'd like to keep QGIS updated. What will happen if I update from 2.8.1 to 3.12.2? I seem to recall updating it in the past and having to put quite a lot of things right manually before I could use it again. Would it be a more straightforward process this time?


texts: Creating Slope Map Using Contours in ArcGIS Desktop?

texts: I am trying to create a slope map. I have many contour elevation files. I tried to use the slope tool in ArcToolbox, but it calls for raster files and I'm working with vector data.


texts: VRTs (vector union layers) not working in QGIS

texts:
I have created a set of VRT's for base mapping vector tile layers following this  tutorial

This works well for displaying the layers but If I want to export or copy out a subset of the features the resulting set of feature is different to what I had selected. When I draw a selection area the selected features don't correspond to the area I have drawn, features 2KM away seem to get selected randomly.

If I do Project> export dxf I get the right features/ area I wanted. There is not an option to create a spatial index.

Maybe vectors in VRT just isn't properly supported.

QGIS 2.10.1 OSX 10.11.1


texts: Corine Land Cover data is not the same as in reality

texts:
I downloaded a CLC2000 shapefile and clipped it after a county shapefile in QGIS.
The goal was to calculate the total land area, but the amount obtained of 983.759 Ha does not coincide with the real area of the county of 446.600 Ha.
I can not figure out the error,if it's from map CRS or from other sources.


texts: Postgis SQL schema for features expected by GeoTools

texts:
I am building a small application which uses postgres and has a couple of tables which have a postgis polygon as one of their columns. I was thinking of using a Geotools DataStore rather than querying it directly through SQL, so that I can take advantage of the library tools (and potentially be able to change to a different DataStore easily if needed.)

I am trying to understand how to get Geotools to know my Schema and which columns it should consider as attributes to my features. All the examples I found in the documentation just set the database URL, username, password, database name etc. but nothing about the tables to extract the features from.

http://docs.geotools.org/stable/userguide/library/jdbc/postgis.html

How do I tell the PostGIS Plugin how to extract the polygons and attributes from my database?

Alternatively, does Geotools expect a table with a certain name and certain structure? If yes where can I find the information about what it expects?


texts: Leaflet Marker placement with minimum distance

texts:
I am using Leaflet and plotting huge marker. There are a lot of places where the markers overlapping each other like this (Don't confuse to see Google map layer. I am using leaflet.)

Is there any option that can show only one marker within a given radius. I am aware of Clustering Plugin. I do not need a cluster. I want to visualize the marker like Tilemill/Mapbox Studio where marker-min-distance plot markers where minimum distance between two markers can be provided. On map zoom again the marker-min-distance will be applied and more marker will be visible.


texts:
I am trying to use an attribute of the data driven page in the legend. I am labeling a polygon layer set to only show if it matches the ddp through a Page Definition.

This is what I want:

This is what I'm trying:

I'm using ArcMap 10.5 with data driven pages enabled.

Has anyone tried doing this before?

I'd like to stay within ArcMap for ease of passing this to someone less comfortable with ArcPy, but would be open to a scripting solution also.


texts: Adding relationshipClass into ArcMap map when cardinality is not M-N?

texts:
I am wondering why it's possible to add a relationshipClass from the catalog into a map in ArcMap 10.3 only when the cardinality is M-N (many-to-many).

When I try to drag and drop other relationshipClass with cardinality (1-1 or 1-M) I have the following error:

  Could not add the specified data object to the map.

Is there any specific reason to that?


texts: How to remove duplicate points based on other attributes

texts:
I have a point dataset of hundreds of thousands of points (zip+4s), many of which share the same value for one attribute (zip+4), but have different values for other attributes (ex: address). In other words, each zip+4 is represented by multiple points at different addresses. I need a dataset that has only one record for each zip+4. In the past, I have simply dissolved the entire dataset by zip+4 and ended up with multi-point features. This, however is starting to cause problems down the line, because I need the location of each zip+4 to be unique - and if a single zip+4 crosses into multiple cities or counties, it gives me issues.

I have gotten as far as dissolving them into smaller multi-point features by city/county, so that each unique combination of zip+4, city, and county is represented by it's own multi-point feature (and I retained a count attribute of how many points make up each multi-point feature). Now how do I get rid of all but the one record for each zip+4 that has the most points in it?


texts: translate point based on the point side attribute and link direction

texts:
I have address points with house number and their side.
For instance, there is a point of house number 10, and it's side is Left (whether Right or Left)

I also have the link that has relationship with the point (it shares same link_id)

what I want to do is translate the location to certain distance away from the link.

I think I probably have to draw perpendicular lines that cross the link at the house address points and move points based on the points side attribute.

I believe the centroid of the perpendicular line is the house address point
the link has from node and to node so I think I can calculate the direction yet I do not know how to find out direction of the link in this type of shape

I do not know where and how to start :(
What I have found so far is
st_azimuth, st_makeline, st_translate..


texts: Converting point sets to polygon boundaries?

texts:
I'm trying to convert sets of points to their respective polygon boundaries.  I think this is something like a Voronoi diagram or Convex Hull, but not quite.  I'm sure there's a technical term for it, but I'm a beginner to GIS.

It's probably best illustrated with the following image:

So, given 4 sets of points (each set in its own colour), what are the four polygons surrounding each respective set of points such that the spacing between polygons maximizes the margins and all the polygons meet together?

Here's a simplified example with 3 points:

A,43.656943,-79.393928,blue
B,43.66663,-79.402682,green
C,43.656447,-79.408004,green

You'll have to view the following graphic for it to make sense:

Essentially, each point has a colour and you expand a buffer of a certain distance out from it to start.

The hard part is when there's a collision with other points (their boundaries) and where to compute that boundary. The result is a set of polygons (one for each colour/set).


texts: EVI not producing values between 1 and -1 Landsat5/7 TM ETM+

texts:
I am trying to do EVI in ArcMap 10.3. I have done the atmospheric corrections for all the bands needed but cannot get the values between 1 and -1. I have used the formula EVI = G *(Band4 - Band3)/(Band4 +C1*Band3 - C2*Band1 +L)
L=1, C1 = 6, C2 = 7.5, and G (gain factor) = 2.5.

Can anyone suggest what the problem might be?

Or if the formula comes in a different format that could work?

Even better, do you have the exact formula to input into ArcMap raster calculater with brackets in the right place?


texts: Question on Exporting data from 'collect values' output in ArcGIS model builder

texts:
I have a question regarding an answer to "Exporting data from 'collect values' output in ArcGIS model builder".  I'm trying to do the same thing as the author to that question and I've followed the steps and the code in that answer.  The code I am using is below:

multivaluesToCsv(r"%CSV File%", "%Feature Class Name%", "%NNZScore Value%", "%NNExpected%", "%NNObserved%", "%PValue%", "%NNRatio%")

import csv

def multivaluesToCsv(csvfile, featureclassnames, nnzscorevalue, nnexpectedvalue, nnobservedvalue, pvaluescore, nnratiovalue):
with open(csvfile, 'wb') as f:
 w = csv.writer(f)
 columnNames = ["FC", "NNZ", "NNE", "NNO", "PV", "NNR"]
 w.writerow(columnNames)
 rows = zip(*map(lambda x: x.split(';'), [featureclassnames, nnzscorevalue, nnexpectedvalue, nnobservedvalue, pvaluescore, nnratiovalue]))
 w.writerows(rows)

My problem is that the CSV file is only being populated with the final row of values.  The titles all appear but only the last iteration's data is included.  Any suggestions?


texts: Finding upstream locations (polygons) using ArcGIS Geometric Network?

texts:
I am New to Geometric Networks.
I need to build a river network and find all locations (polygons) upstream.
This is the template I am trying to use.

I have my rivers with to and from node and the flow direction and accumulation.
I have my pour point (flags) snapped to the network.

How do I build the network to ensure that flow is only downhill and how do I get the upstream IDs?

In the image below I would need to return the Blue Circles (the polygon IDs) for the green pour point. At the moment my network appears to be a free-for-all. When I trace upstream everything is returned.

The network works now. The snapping (for my flags) and the set flow direction worked. My remaining question is how to get the Polygons out! I have 50,000 Flags so I assume I iterate a model, convert the selection to a feature class, then intersect the polygons.

Is this correct or does an easier way to get the polygons out exist?

I feel I need to get the Polygons IDs onto the rivers and then just trace upstream.


texts: Methods for optimizing multcore processing in ArcGIS

texts:
I am interested in learning methods to utilize the full extent of multicore processing power available on a desktop computer.  Arc states that background geoprocessing allows the user to utilize multiple cores, however, tasks essentially have to wait in line for the previous task to be completed.

Has anyone developed parallel or multithreaded geoprocessing methods in Arc/Python?  Are there hardware bottlenecks that prevent multicore processing on individual tasks?

I found an interesting example in Stackoverflow that caught my interest, although it is not a geoprocessing example:

from multiprocessing import Pool
import numpy

numToFactor = 976

def isFactor(x):
    result = None
    div = (numToFactor / x)
    if div*x == numToFactor:
        result = (x,div)
    return result

if __name__ == '__main__':
    pool = Pool(processes=4)
    possibleFactors = range(1,int(numpy.floor(numpy.sqrt(numToFactor)))+1)
    print 'Checking ', possibleFactors
    result = pool.map(isFactor, possibleFactors)
    cleaned = [x for x in result if not x is None]
    print 'Factors are', cleaned


texts:
I am using overpy to collect data information about taxiways, taxilanes and parking positions of an airport (icao code LOWW). I was able to achieve my goal I wrote the following lines of code:

import overpy
api = overpy.Overpass()
r = api.query("""
/*
This has been generated by the overpass-turbo wizard.
The original search was:
“taxiway in airport ”
*/
[out:xml]/*fixed by auto repair*/[timeout:25];
// fetch area “airport” to search in
area[icao~"LOWW"]->.searchArea;
// gather results
(
  // query part for: “aerodrome”
  node["aeroway"="aerodrome"](area.searchArea);
  way["aeroway"="aerodrome"](area.searchArea);
  relation["aeroway"="aerodrome"](area.searchArea);

  // query part for: “taxiway”
  node["aeroway"="taxiway"](area.searchArea);
  way["aeroway"="taxiway"](area.searchArea);
  relation["aeroway"="taxiway"](area.searchArea);

  // query part for: “taxilane”
  node["aeroway"="taxilane"](area.searchArea);
  way["aeroway"="taxilane"](area.searchArea);
  relation["aeroway"="taxilane"](area.searchArea);

  // query part for: “parking_position”
  node["aeroway"="parking_position"](area.searchArea);
  way["aeroway"="parking_position"](area.searchArea);
  relation["aeroway"="parking_position"](area.searchArea);

);
// print results
out meta;/*fixed by auto repair*/
>;
out meta qt;/*fixed by auto repair*/
""")

However my goal is to save the response into the .osm file format. Unfortunately i cannot find in the overpy documentation a way to do it.
Would you be able to suggest a way to achieve my goal?

NOTE: I was thinking to use other packages as well such as osmnx but i cannot get the query as right as with overpy.

import osmnx as ox
G = ox.graph_from_place('wien flughafen austria', infrastructure='way["aeroway"]',simplify=False, custom_filter = ('["aeroway"~"taxiway|taxilane|parking_position"]'))
ox.save_load.save_as_osm(G, filename='filename.osm')


texts: Finding closest point to shapefile coastline Python

texts:
I have a dataset consisting of 2,2 million houses in a Norway, and I want to find the distance to the coast for each dwelling. The location of the house is made up by latitude and longitude coordinates, and I have downloaded a shapefile of Europe with all the coastlines from (Could not find a shapefile for only Norway):
https://www.eea.europa.eu/data-and-maps/data/eea-coastline-for-analysis-1/gis-data/europe-coastline-shapefile

Which have the coastline in both MULTILINESTRING and POLYGON format.

Then I try using the nearest point function from the shapely library, but it only assign the same value for each point, and it takes too long to run.

import geopandas as gpd
from shapely.geometry import Point
from shapely.ops import nearest_points

def findClosestCoastline(df):

   #assign Point from the coordinates of dwelling
   df['point'] = [Point(x, y) for x, y in zip(df['coord_x'].astype(float),
                                              df['coord_y'].astype(float))]
   #Load shapefile with MULTILINESTRING format
   map_2= gpd.read_file("Europe_coastline.shp")

   for i in range(df)
      #TODO: Fix this, get the closest point from the MULTISTRINGLINE
      np = nearest_points(map_2.iloc[0]['geometry'], df.iloc[i]['point'])[0]
      df.loc[i, 'closestPoint_Coast'] = np

      #TODO generate the distance in KM
      df.loc[i, 'distance_To_Coast'] =

return df


texts: Calculate Geometry function returns underestimated distances

texts:
I'm making a map of the Colorado Trail and obtained a set of waypoints taken down the entire trail. I projected the set of waypoints from WGS84 (unprojected) to NAD 83 UTM13N after adding to the data frame which is also in the latter projection). Coordinates display as  northing/easting meters as expected. Everything looks fine but using Calculate Geometry to compute trail length (total and individual section lengths) returns erroneous underestimates of published lengths (the published lengths have been measured with a distance wheel). Not even close-- off by nearly 100 miles for the whole trail.
I computed surface length using a 1/3 arc sec DEM but this in no way accounts for the error. I can't figure this out but I feel it has something to do with the projection and how it relates to the initially unprojected waypoints.
Can anyone clear this up as I am out of ideas?


texts: ArcGIS Javascript API returning value to popup from function

texts:
I'm creating a custom pop-up where I send the attribute table value into a function to reformat it. It's called like so:

var herbicide = '${Herbicide1:herbName}'

It then enters the function like this (it has many more cases than what's shown):

herbName = function(value){
   herb = '';
   switch(value){
       case null:
           herb = 'None';
           break;
       case 1:
           herb = '2,4-D';
           break;
       ....
   };
   return herb;
};

The 'return herb' line should return the new value and I thought that the 'herbicide' variable would be set equal to whatever was returned, but that does not appear to be the case if I actually try to print it's value. When I print it:

console.log(herbicide);

It returns
${Herbicide1:herbName}

This is making it really hard to do conditional statements on the variable. For instance, if herbicide is set to 'None' then I don't want it to show up in the popup, but using this statement:

if (herbicide == 'None'){
    //do something
}

It will never evaluate true because unfortunately, no matter what value was returned from the function, the value being tested against is always the string '${Herbicide1:herbName}'

What's weird to me is that you can't see the value by printing it to the console, yet I know my functions are returning the correct values because it prints them correctly to a popup. I feel like I'm doing something stupid or backwards, anybody know what I need to do?

Using JS API 3.x


texts:
When a leaflet popup shows up, I would like to animate its offset value, so that the popup moves off to the right when it pops up. I imagine the codes goes something like this, but I haven't been able to figure it out:

layer.bindPopup(content,
{
offset: [x, y]
});

layer.on('click', function () {
setInterval(functionThatAnimatesxandy, 1000);
}

So apparently you can not animate the offset value. What might be a solution with CSS transitions?


texts: Exact conversion of ECEF to geodetic coordinates

texts:
I came across the paper "Exact conversion of Earth-centered Earth-fixed coordinates to geodetic coordinates", which describes a series of equations which can be used to calculate a set of geodetic coordinates from a set of ECEF coordinates. However, there seem to be a plethora of methods for doing this type of calculation, and so I'm wondering:

Is this method mathematically correct?
Does there exist a closed-form solution to the conversion problem which is widely accepted as a standard in the GIS community?


texts: Smoothing/Aggregating polygon in ArcGIS Desktop?

texts:
I have this one big polygon separated by a big gap/space like this

Can anyone advise how can i edit, so that the space should be narrow, smooth without a big gap like this?

What function should i use? cartography tools? Spatial Analyst function? Watershed?


texts: UpdateCursor not filling correctly with conditional statement

texts:
I am attempting to update a table so that a field, "ALLREADY", fills with "Y" or "N" depending upon if "dmready" equals "total_records" AND if "dm_stat" is equal to "Ready for DM".  When I run my script, all records fill with N, even those that meet the conditions (as you can see from the first record in the screenshot below).

My script:

inTable = os.path.join(env.workspace, "pivotTableTEMP")
with arcpy.da.UpdateCursor(inTable, ["dmready", "total_records", "dm_stat", "ALLREADY"]) as cursor:
    for row in cursor:
        if row[0] == row[1] and row[2] == "Ready for DM":
            row[3] = "Y"
            cursor.updateRow(row)
        else:
            row[3] = "N"
            cursor.updateRow(row)

I have double checked that all fields in my script match the column names and that Ready for DM does not have any hidden characters.


texts: How to generalize and vectorize a LandUseClassification raster file?

texts:
I want to vectorize that rasterfile with its 6 landuse classes. The classification is very noisy and have to be generalized before.

The result should be a polygone vector file.

I use SAGA, (GRASS) and QGIS.


texts: Shapefile export of holes geometries in QGIS with wrong type MultiPolygon

texts:
I want to export a layer with holed geometries in it. But when I save my features to shapefile in QGIS 3 it creates a file with MultiPolygon as geometry type. The fact is that a holed polygon isn't a multi-surfacic polygon and I want to export it as Polygon only.

I tried to modify manually the WKT with geom_from_wkt() with following formula :

geom_from_wkt( replace( replace( geom_to_wkt( $geometry ), 'MultiPolygon (', 'Polygon '), ')))', '))' ) )

With Geometry by expression and the temporary layer created by this algorithm is correct (with Polygon type) but when I save it to shapefile it puts MultiPolygon back again.


texts: Using ASCII LiDAR data tiles which are centimeter rounded (n Z) in ArcGIS Desktop seems to round them to meters?

texts:
I try to use four LIDAR tiles in .asc with arcgis9.3. All of this tiles are centimeter round (in z) but one of them seems to be round to meter (see tile on the left up corner on the image below), while opening it with notepad we can see that it is centimeter round.

I compared the header of that four files and i didn't find any difference.

On that dropbox link you can find two LIDAR tiles. One is correctly read by arcgis when the second is not.

https://www.dropbox.com/sh/1pwuqr0nppb42qu/AAAk0zeVHbpVO1j9KZKMBf7Ca?dl=0


texts: python rioxarray.clip masking netcdf data with a polygon returns all nan

texts:
My question is similar to what others have already asked but the posted solutions haven't worked for me. The similar posts are masking a netcdf file using a shapefile of points with rioxarray and how to mask netcdf time series data from a shapefile in python.
I am simply trying to clip an xarray DataArray with a polygon using rioxarray. It seems the other posted solutions should be working for me but when I clip I either get all nans or no clipping at all. Can anyone see what I am doing wrong? Hoping this is just a derp mistake somewhere
import xarray as xr
import geopandas as gpd
import rioxarray as rio
from shapely.geometry import mapping, Polygon
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# create contiguous US polygon for clipping CMIP5 data
world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
US=world[world['name']=='United States of America']  # conus, AK, HI, PR
conus_poly = Polygon([(-126, 24), (-126, 50), (-66, 50), (-66, 24), (-126, 24)]) # conus box
conus=gpd.clip(US,conus_poly) #cut out AK, HI, PR
conus.geometry

I have also tried reading in a conus shapefile instead of the above method but get the same results.
# load CMIP5 data and add crs info
lat1,lat2,lon1,lon2=15.,55.,230.,300.
modDS=xr.open_dataset('/mnt/d/DATA/CMIP5/Amon/pr/regridded/regrid_pr_Amon_ACCESS1-0_historical_r1i1p1_190601-200512.nc').sel(lat=slice(lat1,lat2),lon=slice(lon1,lon2))
pr=modDS.pr
pr.rio.write_crs("epsg:4326", inplace=True)

The data file I'm using above can be accessed here
# clip and plot
clipped1 = pr.rio.clip(conus.geometry.apply(mapping), conus.crs, drop=False, invert=False)     # clips everything - all nan
clipped2 = pr.rio.clip(conus.geometry.apply(mapping), conus.crs, drop=False, invert=True)      # clips nothing

plottime='2000-08-01'
fig=plt.figure(figsize=(10,10))
ax1=fig.add_subplot(211,projection=ccrs.PlateCarree())
clipped1.sel(time=plottime).plot(ax=ax1, levels=[0,1,2,3,4,5,6,7,8,9,10])
conus.boundary.plot(ax=ax1,color='black')
plt.title("clipped pr invert=False")

ax2=fig.add_subplot(212,projection=ccrs.PlateCarree())
clipped2.sel(time=plottime).plot(ax=ax2, levels=[0,1,2,3,4,5,6,7,8,9,10])
conus.boundary.plot(ax=ax2,color='white')
plt.title("clipped pr invert=True")

plt.show()


texts: Customized PDF Report of Each Manhole from a Database

texts: I need a way to generate individual PDF reports for all the data collected for each manhole in a database.  I would end up with hundreds of reports, one for each manhole in the feature class.  The report needs to selected attribute data such as manhole type in a chart along with 1-3 pictures and possibly a small map.  It doesn't look like the reporting tool built into ArcMap can make something this complicated.  Any help would be greatly appreciated!


texts:
I built an extension module for Geoserver 2.5 using Maven and Wicket. It places a Label to the left of the navigation list that calls a HTML-Page in Geoserver when clicked. It looks and works as it should; the problem is that it is visible and usable without login.

How can I configure my extension module/GeoServer so that it gets to the plugin only after login?


texts: Bolding Text String Not Working in ArcMap

texts:
I am using ArcMap 10.5 with an Advanced level license.
I am wanting a portion of my text box to be bold. I clicked Insert --> Text. When I type <bol> it literally shows up on my text box as <bol>. I have searched everywhere for recommendations, tried adding quotation marks and nothing works.


texts: How to refresh the view in QGIS Print Composer using Python?

texts:
I am trying to automate the printing of maps based on different attribute classifications using Python in QGIS. Everything is okay except for when I open the printed image, the labels are blurry. Manually, I can fix this problem by clicking the "Refresh View" button in the print composer window, but I cannot find the right way to automatically refresh the map in the print composer based on the QGIS API. Here is printing part of this code:

def printMyMap(canvas, field, itemp):
    canvas.refresh()
    mapRenderer = canvas.mapRenderer()
    c = QgsComposition(mapRenderer)
    myFile = 'C:/Users/User1/'+itemp
    myTemplateFile = file(myFile, 'rt')
    myTemplateContent = myTemplateFile.read()
    myTemplateFile.close()
    myDocument = QDomDocument()
    myDocument.setContent(myTemplateContent)
    c.loadFromTemplate(myDocument)

    #Map
    cmaps =  c.composerMapItems()
    cmaps[0].updateCachedImage()
    cmaps[1].updateCachedImage()

    #Legend Titles
    #legend = c.getComposerItemById("3")
    #legend.setTitle(ltitle)

    #Show Composer
    #iface.actionPrintComposer().trigger()
    composerList=iface.activeComposers()
    composerView=composerList[composerList.index(max(composerList))]
    composer = composerView.composerWindow()
    composer.show()

    #Save out
    dpi = c.printResolution()
    dpmm = dpi / 25.4
    width = int(dpmm * c.paperWidth())
    height = int(dpmm * c.paperHeight())

    image = QImage(QSize(width, height), QImage.Format_ARGB32)
    image.setDotsPerMeterX(dpmm * 1000)
    image.setDotsPerMeterY(dpmm * 1000)
    image.fill(0)

    imagePainter = QPainter(image)
    sourceArea = QRectF(0, 0, c.paperWidth(), c.paperHeight())
    targetArea = QRectF(0, 0, width, height)
    c.render(imagePainter, targetArea, sourceArea)
    imagePainter.end()
    image.save("F:/Maps/"+field+".png", "png")


texts: Connecting QGIS to ArcGIS world imagery?

texts:
Have connected QGIS to ArcGIS World Imagery using python console as described on http://www.northrivergeographic.com/qgis-adding-an-arcgis-rest-service.  QGIS Seems to connect and make a raster layer the image is not correct.

I'm very new to QGIS.


texts: Is there a way to perform the same task as the gdalbuildvrt utility using GDAL Python bindings?  So far I have not found any way to do this other than creating a vrt of a single dataset and manually editing the xml.  I would like to create a vrt from multiple rasters (essentially performing a mosaic).  Is this possible using pure Python?  My other option is to use subprocess to simply call gdalbuildvrt.


texts:
I am using GDAL to convert GeoTIFF to KMZ. I am using the following command:

gdal_translate.exe -of KMLSUPEROVERLAY *.tif out.kmz

If I have one TIF file in the folder, it works ok. If I have more than one, then I receive the following error:

FAILURE: Too many command options.

How can I run the gdal_translate.exe command to combine more than one TIF file into one KMZ file?


texts: Missing label points with fishnet creation

texts: I am trying to create a fish net with label points in Arcmap 10. I can get the polyline or polygon of a grid on my project, but I am missing the label points on the map in both data and layout view. In creating the fishnet I have clicked the create label points box. A layer is created and shows up in my TOC, but when I zoom to layer it takes me to white space. The attribute table for the point label layer is also empty.


texts: Write GeoJSON As Feature Object in Postgis

texts:
I'm trying to create a geofence/circle from a lat/long set and I'm close except that the GeoJSON that PostGIS creates differs in structure from the client library that our front-end uses.  Here is my python to create the geofence:

session = DatabaseSessionService().create_session()

METERS_IN_HALF_MILE = 805

longitude = 21.304116745663165
latitude = 38.68607570952619

geom = func.ST_Buffer(func.ST_MakePoint(longitude, latitude), METERS_IN_HALF_MILE)
geofence = Geofence(
    geom=geom,
    geojson=cast(func.ST_AsGeoJSON(geom), JSON),
    company_id=188,
    name='Test Fence 1'
)

session.add(geofence)
session.commit()

This will create a record with geojson set to:

{
   "type":"Polygon",
   "coordinates":[
      [
         [
            826.304116745663,
            38.6860757095262
         ],
      ...
      ]
   ]
}

However I'd like to maintain the convention set previously and have the GeoJSON written as a Feature Object like this:

{
   "type":"Feature",
   "properties":{

   },
   "geometry":{
      "type":"Polygon",
      "coordinates":[
         [
            [
               -92.29188859462738,
               38.913722221672636
            ],
            ...
         ]
      ]
   }
}

I've looked into using ogr2ogr but I'd prefer that this be an automated process that happens on creation as opposed to being a post-creation task.


texts: Click event doesn't work with Google Maps

texts:
I'm trying to create a simple click event when I click on a list item on my page, a marker should be created on the google map in that page. But it does not happen. Is there some google maps API limitation I'm missing here?

<!DOCTYPE html>

  <meta charset="utf-8">
  <title>Google Maps Test</title>

Roman Weekend

Termini

and here is the script:

function initMap() {
  // Create a map object and specify the DOM element for display.
  var map = new google.maps.Map(document.getElementById('map'), {
    center: {lat: 41.902728, lng: 12.496246},
    scrollwheel: false,
    zoom: 13
  });
}

google.maps.event.addDomListener(window, 'load', initMap)

  var marker = new google.maps.Marker({
    position: {lat: 41.900819, lng: 12.501757},
    title:"Termini"
});

function terminiFunction() {
  marker.setMap(map)
}


texts: OpenLayers with Yandex layer position issue

texts:
I am kind of new for openlayers and I am working on a project that have many different layers.
one of this layer is Yandex, I added this layer few days ago and it seems to work fine, untill QA :) discovered that when you switch from this layer to another it not position the same as the other layer, I added this layer base on the following article (couldn't find anything else that show how to add the Yandex)
http://xbb.uz/openlayers/i-Yandex.Maps

when I have lines layer on the map and change the base layer from google, bing or any other layer to yandex the lines or the yandex layer are not positioned right.
YANDEX layer

OPEN STREET MAP LAYER


texts: How can I create a polygon(filled polygon) from a single pixel stored in a hdf file(MODIS)?

texts:
I'm working with MODIS data to get EVI. I'm following this question to get the required values of EVI. The hdf file has 4800x4800 rows/columns (250m). I can get the X,Y pixel coordinates and the corresponding value. Now I'm trying to get this pixel as a polygon so that I can merge all the polygons with the same value. But I'm not sure how to get it.

As per my understanding, the X,Y pixel values are the coordinates of bottom right corner of the pixel. How can I create a polygon out of it ?


texts: ModelBuilder: How to extract data by date, export to separate tables and interpolate tables to rasters

texts:
I've been using ArcGIS for a couple years, but when it comes to model builder I've only worked with very simple models.  I'm hoping someone that is much smarter than myself can help me in building a complicated model.

I've got a very large table of data (about 300k rows) that contains points with lat/lon data, wind speed and date/time.  There are about 60 different points with wind values that have been recorded for each hour during 2012.  What I'd like to do is build a model that will:

Iterate through each row and select/export the data associated with first hour (00Z) Jan 1st 2012 (about 60 rows) to a separate table, then the 2nd hour (01Z) of Jan 1st 2012 to its own separate table and so on through all of 2012.
Perform Inverse Distance Weighting on all of the hourly tables with lat/lon and wind values that were created from step 1 above
Reclassify the rasters so that points in the raster that have wind values of 5 or greater are reclassified to "1" and all other cell's wind values are reclassified to "0"
Then iterate through the rasters created in step 3 and use something like "Collect Values" & "Cell Statistics" to sum up all of the rasters into a final summed raster for all hours in 2012.

Was going to include screenshot of my data, but I guess since this is my first post I'm not trusted yet.  I've broken the date down several different ways (just date, just time, date-time, and date-time in text format) because I've had issues trying to iterate through and select the data in my original, very large table only by date.

Hopefully I provided enough information, but if not I will try to go into more detail.  Thanks in advance for any assistance.


texts: Why won't my hyperlinks work in the CartoDB infowindow?

texts:
I've been setting up a map of project locations for my company, and I want the infowindow to include a link to the more detailed project description page on our website. I have checked several of the other questions related to URLs in CartoDB infowindows, but none seem to cover the same issue.

I uploaded a dataset with the URLs for our project pages in a field called "Link" and I tried using the following in the Custom HTML window (along with a thumbnail image of the project and a couple other details):

more info

When I check the "more info" link in the infowindow on the map, it always points back to the map itself (https://bfradkin.cartodb.com/viz/ce4b4070-87f8-11e5-b01c-0ea31932ec1d/map), not the URL stored in the Link field.

See below for a screenshot showing the map URL when I hover over the link:

I'm hoping to solve this without getting into CartoDB.js, any thoughts?


texts: I'm trying to change my coordinate formats in RStudio

texts:
So my data entries for Longitude and Latitude look like this:
87010454" "85404107" "86294685" "85545624" "86243217" "86135700"
And I am trying to plot a map of the US and add my data to it using this codes:

library(ggmap)
library(RCurl)
library(xlsx)
urlfile<-'http://www.psc.isr.umich.edu/dis/census/Features/tract2zip/MedianZIP-3.xlsx'
destfile<-"census20062010.xlsx"
download.file(urlfile,destfile,mode="wb")
map<-get_map(location='united states', zoom=4, maptype="terrain",source = 'google',color = 'color')
gmap(map) + geom_point( aes(x=LONGITUD, y=LATITUDE, show_guide = TRUE, colour=STATE), data=accident, alpha=.5, na.rm = T) + scale_color_gradient(low="beige", high="blue")

However I think the format for this variables is wrong because I am able to get a map, the data won't appear on it though. Is there any way I can change my format in R Studio to get it to work?


texts:
I teach a non-geographical technical course, and the topic will be triangulation (..of suspects via mobile phone signals while using an archived well-known case). Since the system in questioned is GSM, we/they used the TA value (distance from tower quantized from [0:63], where each step is ~550 meters). So if TA is (eg.) 10, the distance from the tower is (in an ideal case) from "10*550m" to "10*550m+549.99m".

I have coordinates of three towers in lat,lon (WGS84), I have the TA values, and I would like to draw three donut shapes (inner radius=ta*550m, outer (ta+1)*550m) which would ideally intersect in one relatively small area. I also have a raster map (satellite images). I have already verified the raster layer and the basestation points are visualised OK, but I have no idea how to draw the 'donuts' (except photoshop and eyeballing it).

What would be the simplest way to visualise that in QGIS?


texts: In ARCGIS 10.2 can I create a 2 page map document that prints an image on the back side of the page?

texts: From arcgis 10.2 I need to create a pdf file that will print a map on the front side of a page and an image on a back side of the page. How would I do this?


texts: Using Leave one out method on kriged data

texts:
I am looking for a function to do the following.

Given a boundary file, and soil sample points within the boundary I am wanting to use the Leave one out validation technique.

What I mean is create a Kriged map leaving out n point, and then compare the kriged value with the known value.

The purpose of this is to identify areas within the boundary that need resampling if the kriged value is significantly different from the actual value.

Does anyone know of a method within QGIS or custom python script to do this?


texts: Convert ArcGIS SWM to numpy array directly

texts:
How do I convert an ArcGIS (10.3) spatial weight matrix (SWM) file to a numpy array directly?

swm -->(swm to table)--> table -->(table to numpyarray)--> numpyarray

is very slow with 16.000.000 entries. There must be a way to directly (or at least faster) convert a SWM to numpy array, is there not?


texts:
I'm a new users of spatial analyses packages in R.
I have created a SpatialPolygonsDataFrame which contains a number of rectangular polygons that i named "Square1" to "Square64".
Until now there is no problem i can plot these polygons on a map. The problem is that i have a data.frame that contains the number of vesels in each "Square_i" (i=1 to 64) polygon, but i don't know how to relate each polygons in my SpatialPolygonsDataFrame object to the data.

Can anyone help me?


texts: bulk save mapinfo tab files with new projection

texts:
Is there a way to builk convert mapinfo files to a new projection (and output in the same mapinfo format) without using the mapinfo sdk?

I did look into ogr2ogr but I couldn't get it to work even when I added the to/from epsgs that I needed to the csv files

I need to convert from EGSP 28356 -> EPSG 3857


texts: LiDAR Data captured, unsure what units or format this time field is in?

texts: I haven't been able to figure out or find out from searching around what units or format this time field is in my LiDAR data that was collected?


texts:
I want to load a GeoTIFF file from a browser using Leaflet. I have created a custom control button that on click loads a file from local computer.

fileInput.Change() is the function that is called, when a user select a tiff file. Right now I am loading a file from local directory, and it's size is around 57 MB, and I am getting error of 206 in my browser.

    // load raster button
    var customControl =  L.Control.extend({

      options: {
        position: 'topleft'
      },

      onAdd: function (map) {
        var container = L.DomUtil.create('input');
        container.type="button";
        container.title="No cat";
        // container.value = "42";

        container.style.backgroundColor = 'white';
        //container.style.backgroundImage = "url(https://t1.gstatic.com/images?q=tbn:ANd9GcR6FCUMW5bPn8C4PbKak2BJQQsmC-K9-mbYBeFZm1ZM2w2GRy40Ew)";
        container.style.backgroundSize = "30px 30px";
        container.style.width = '30px';
        container.style.height = '30px';

        // container.onmouseover = function(){
        //   container.style.backgroundColor = 'pink';
        // }
        // container.onmouseout = function(){
        //   container.style.backgroundColor = 'white';
        // }

        container.onclick = function(){
        document.getElementById('fileInput').click();
          console.log('buttonClicked');
        }

        return container;
      }
    });

    fileInput.onchange = function(e) {
      debugger;
      var bounds = new L.LatLngBounds(
            new L.LatLng(-179.9999999999999716,-89.9999999999820091 ),
            new L.LatLng(179.9999999999640465, 89.9999999999820091)
          );

            // var overlay = new L.imageOverlay("GRAY_50M_SR/GRAY_50M_SR.tif", bounds)
            //   .addTo(map);
  //           var  bounds = new L.LatLngBounds(
        // new L.LatLng(40.71222,-74.22655),
        // new L.LatLng(40.77394,-74.12544));

            var overlay = new L.ImageOverlay("GRAY_50M_SR/GRAY_50M_SR.tif", bounds, {
                        opacity: 0.5,
                        interactive: true,
                        attribution: '&copy; A.B.B Corp.'
                    });
                    map.addLayer(overlay);

    };

How can I load a GeoTIFF file in
Leaflet from browser?


texts:
I have a table with ~15K records consisting in points with a 200m buffer and I need to investigate how many of them are in proximity of polygons contained in another table (~20K).

I wrote a function (really simple one) which uses a record variable (rowgeo) and then uses ST_Overlaps like this:

SELECT gid
FROM pi a
WHERE ST_Overlaps(a.geom, rowgeo.buffer) INTO var

well, the function is dramatically slow... it's running since yesterday morning. Both tables have a gist index on geometry column.

Any hints?


texts: Multiple layers auto-updated when one of them is modified

texts:
I just discovered this post and I found what I searched here. But I was wondering if there was a way to have multiple layers (let's say l1, l2, l3, l4) being updated automatically when you change something on one of them.

For example if I change the outline style of l1 from solid to dot I would like to see l2, l3 and l4 updated the same way without copy/pasting the style like mentioned in the link

(I'm using qgis 2.18.28 on windows)


texts: How to fit map view to a layer extent in OpenLayers5?

texts:
I added a number of markers to the map:

const map = new ol.Map({
    target: 'map',
    layers: [
        new Layer.Tile({source: new OSM()}),
        markersVector,
    ],
    view: new ol.View({
        center: [0, 0],
        zoom: 2,
    })
});

// ...

let setMarkers = function (markers) {
    let features = markers.map((marker) => createFeature(...marker));
    markerSource.addFeatures(features);

    // how to fit?
};

setMarkers(window.markers || []);

After adding the features I would like the map to zoom as much as possible so that all markers are visible at once. How to do that? Previous versions had a zoomToExtent method which is now missing, obviously.

I tried the map.getView().fit(...) method with various inputs, but I didn't manage. The view stays fixed to such location:

I tried, but it does not do anything:

map.getView().fit(
    map.getView().calculateExtent(),
    {constrainResolution: false, padding: [10, 10, 10, 10]}
);

Edit:

The map.getView().calculateExtent() returns the same extent initially and after adding the features. So I guess my question is how to calculate the extent properly.


texts: QGIS - Calculate net area of classes after raster classification

texts:
I searched on the net but could not get the specific answer. Request for the same.

I have slope map of the area (in percent, derived from r.slope.aspect). I need 4 classes, as, area of slopes less than 5%, 5 to 15, 15 to 20 & above 20%. I got it done using singleband psuedocolour. Now I need the net area for these 4 classes.

How to achieve this in Qgis? Raster to Vector conversion should be avoided as possible. r.stats or r.report gives huge reports. It counts the individual cell instead of class & process accordingly, I think. I only need area for above 4 categories.

Or should I generate all together new raster maps, using raster calculator and with expression as Rasterlayer@1 < 5, etc... & again merge/ integrate into one final..


texts:
I am trying to make a QGIS plugin using Python.
Right now I need way add an image to canvas but I don't want that image to be a new layer in qgis project. Is there a way to do it?

Now I have:

self.canvas = qgis.gui.QgsMapCanvas()
self.canvas.show()

But it shows only white rectangle.


texts: I use save as methods as well as MXD Doctor but nothing happening. I have previously 10.5 version and I have saved my work in MXD. Now I have 10.3.2 version. How can I open .MXD files?


texts: Is this Fastest, most Elegant and Extendable to nD implementation of AABB Intersection?

texts:
I recently did an AABB (Axes Aligned Bounding Box) Intersection implementation on 2D and 3D objects for which I came up with such a short (and very fast) solution as follows.

The demonstration below is for 2D rectangles, nevertheless the algorithm promises extension to 3D, and nD spaces as well, with no additional effort.

AABB has been a standard part of many small and huge computational geometry packages in order to speed up intersection and collision detection, for years. For me, not fan of large complicated libraries, I preferred however to do it myself with no consultation with any of them.

With the above introduction, it is quite possible that I might have missed some situations on those the following simplistic algorithm might fail.
If there was no such failure this algorithm looks like the fastest intersection test for AABB and the most easy-to-extend to nD as well as it requires only few simple numeric comparison.

Introduction:

An AABB is usually defined by two points as (mins,maxs), e.g., in 2D (0,0,1,1) where (0,0) are the minimums for X and Y axes, and (1,1) are the maximums for X and Y axes respectively. Similarly for nD it could be as follows.

P = (x1,y1,z1,...,x2,y2,z2,...)

Let's P and Q be two AABB:

P(mins,maxs)                  % e.g., in 2D (0,0,1,1)
Q(mins,maxs)                  % e.g., in 2D (0.1,-1,3,2)

Intersection between the two P and Q happens only if:

P(mins) <= Q(maxs) & Q(mins) <= P(maxs)

So for the 2D example given above,
P(mins) = (0,0)  and P(maxs) = (1,1)
and for Q:
Q(mins) = (0.1,-1)  and Q(maxs) = (3,2)
Now if we check the conditions mentioned above:
(0,0) <= (3,2) i.e., (0 <= 3) and (0 <= 2) and
(0.1,-1) <= (1,1) i.e., (0.1 <= 1) and (-1 <= 1), we see that all conditions satisfy which means the two AABB P and Q intersect.

Therefore the algorithm requires only few simple comparisons which can be implemented in various forms efficiently. Furthermore, the algorithm can easily be extended into nD.

I wonder if this such simple algorithm is the fastest, the most elegant and extendable to nD for finding intersection between AABBs?

Any comment about potential shortcoming or any discussion or share of experiences would be helpful.


texts: OSM Building as Multipolygon not rendering as building

texts:
I have been working to improve OpenStreetMap data in my local area and I've found something that I think is a problem but I can't figure out how to fix it.

The basic problem is that the Jefferson County Courthouse building in Golden, Colorado doesn't render the building outline on OpenStreetMaps.org, see:  https://www.openstreetmap.org/#map=17/39.72769/-105.20104&layers=N

When I download that bbox in JOSM however, the building outline shows up as a multipolygon named "Jefferson County Courts", with only 1 single part to the multipolygon.  I've created a note in OSM about this, see: https://www.openstreetmap.org/note/428842#map=18/39.72784/-105.20062&layers=N

First question:  Is there a good reason for that building to be a multipolygon?  It is one single building, so I can't think of one.

Second question:  Is there a way to preserve the history associated with the multipolygon, but to just make it a regular polygon?  I'd like to keep the current boundary and the history associated with it, a couple people seem to have put some good time into that shape and I don't want that to be lost.

If converting it to a singular polygon isn't appropriate, what is keeping it from rendering in the default view of OSM?  I know we're not supposed to tag for the renderer, but it's tagged as building=public but isn't showing up like a building.


texts: Applying styling to KML (from a shp) using ogr2ogr?

texts:
Is it possible to generate a styled KML using ogr2ogr (ie pen width and colour?) I cannot seem to find anything in the ogr2ogr documentation.

My code so far is;

ogr2ogr.exe -f "KML" prow-kml.kml "C:\prow.shp"

Unfortunately this just gives me a default colour line when viewing in Google Earth.


texts: Debugging RuntimeError: A column was specified that does not exist from ArcPy?

texts:
My goal is to count the number of wells in each county. I've verified the columns are correct as listed in the attribute table. Do I have a typing or formatting error? I wanted to pull county names from a list but I'm trying to T/S with one variable first.

import arcpy
arcpy.overWriteOutput = True
arcpy.env.workspace = "E:\\GIS4080\\Week5\\Lesson5_Data"

# Variables
wells = "Wells.shp"
counties = "COUNTIES.shp"
Wells_Intersect = "Wells_Intersect.shp"
Wells_Intersect_Layer = "Wells_Intersect_Layer"
input_f = ["Wells.shp", "COUNTIES.shp"]
Date_list = []
wellCnt = 0

# Process: Intersect
arcpy.Intersect_analysis(input_f, Wells_Intersect, "ALL", "", "INPUT")

# Process: Make Feature Layer
arcpy.MakeFeatureLayer_management(Wells_Intersect, Wells_Intersect_Layer)

# List County Names
countyList = [row[0] for row in arcpy.da.SearchCursor(counties, "COUNTY")]
#print countyList

# Count Wells in each County
#for cname in countyList:
cname = "YUMA"
with arcpy.da.SearchCursor(Wells_Intersect_Layer, "COUNTY", '"COUNTY "= ' + cname) as cursor:
    for row in cursor:
        wellCnt = wellCnt + 1
        print cname, str(wellCnt)

  Traceback (most recent call last):
   File "E:\GIS4080\Week5\L5_hw.py", line 32, in
    for row in cursor:
  RuntimeError: A column was specified that does not exist.


texts: Intersecting lines to get points in geopandas

texts:
I'm trying to use GeoPandas to generate points of intersection between two linear shapefiles. For example, I have a shapefile of pretend pipelines (red) and some crossings (blue).

What I would expect, are points wherever they intersect (green):

However, what I'm getting is just two points like so:

Here's my code:

import geopandas as gpd
pipelines = gpd.read_file(r"C:\Temp\GeometryTesting\Pipeline.shp")
crossings = gpd.read_file("C:\Temp\GeometryTesting\Linear_Crossings.shp")

I would expect this to be the correct code to calculate crossing points:

pipelines.intersection(crossings)

or even this:

crossings.intersection(pipelines)

However, both of them return the same result:

0     POINT (147.010039196221 -26.40843953078755)
1    POINT (150.4550260127808 -25.61354916104867)
2                                              ()
3                                              ()
dtype: object

Why am I only getting two points (and two records with no data)? Using geopandas, how can I generate all the crossing points?


texts: ArcMap 10.1 Intersect Not Adding Points at All Intersections

texts: I have run the Intersect tool on a streets centerline layer.  I tried it on a dissolved version of the lines and also a split version of the lines.  The tool is not giving me a point at each intersection in either case.  I have included an image below.  The red circles are an example of the areas that should also have points.


texts:
I am looking for a good parser for the latest GML (Geographic Markup Language) version in preferably: C++/Java/Ruby/...

I am especially interested in reading the geometry of individual features. Does anybody have recommendations?


texts: Vincenty or (Project + Euclidean) for accurate distances

texts:
I posted a similar question on stackoverflow, however was directed here for more help.

If my main goal was to calculate distances such that my error was minimised to the nearest metre assuming that the points I am working are bounded on a country-level (for example UK or France ... those kind of medium-sized countries).

Assuming I start with geographic WGS84 co-ordinates would it be better to:

1) Use the vincenty-formula on the epgs:4326 coordinates

2) Project the co-ordinates (e.g. from epsg:4326 to epsg:27700 for the UK) and calculate the hypotenuse

The latter option would be much quicker; however I am not sure which one would be more accurate.

I ran a nearest-neighbour test using (both unvectorised, haversine rather than the more accurate vincenty ... ideally I would compare both functions vectorised using numpy):

nplen = 1000000
# WGS84 lat/long
point = [51.349,-0.19]

# This contains WGS84 lat/long
points = np.ndarray.tolist(np.column_stack(
        [np.round(np.random.randn(nplen)+51,5),
         np.round(np.random.randn(nplen),5)]))

def proj_list(points,
              inproj = Proj(init='epsg:4326'),
              outproj = Proj(init='epsg:27700')):
    """ Projected geo coordinates"""
    return [list(transform(inproj,outproj,x,y)) for y,x in points]

proj_points = proj_list(points)
proj_point = proj_list([point])[0]

# Run tests
from math import *
def haversine(origin,
              destination):
    """
    Find distance between a pair of lat/lng coordinates
    """
    lat1, lon1, lat2, lon2 = map(radians, [origin[0],origin[1],destination[0],destination[1]])
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
    c = 2 * asin(sqrt(a))
    r = 6371000  # Metres
    return (c * r)

def closest_math_unproj(points,point):
    """ Haversine on unprojected """
    return (min((haversine(point,pt),pt[0],pt[1]) for pt in points))

def closest_math_proj(points,point):
    """ Simple angle since projected"""
    return (min((hypot(x2-point[1],y2-point[0]),y2,x2) for y2,x2 in points))

closest_math_unproj - 22.1 seconds with a distance of 138.491

closest_math_proj - 3.66 seconds with a distnace of 138.904

Where the two points were:

(51.34892,-0.18801)

(51.349,-0.19)


texts: How to split line at a specified distance and getting coordinates of start and end points of each line segment？

texts:
I tried v.split in GRASS：

and my notebook ran a lot of time but showed this error：bad allocation, then I click "OK", and got nothing.

Are there other ways to split lines？and I hope to get 2 new columns in the same shapefile about coordinates of start and end points of each line segment.


texts: Dissolve FeatureCollection using attribute

texts:
I have a feature collection with a number of properties, all numbers. I want to aggregate all features in the collection that have a property in common. For instance, I want to merge features that have the same date value, and I wish to sum up the values of all other fields in the resulting feature collection.
For example, a dummy feature collection:
var feat_list = [
  ee.Feature(ee.Geometry.Point(77.58, 13), {date:'2020-10-02', count:1}),
  ee.Feature(ee.Geometry.Point(77.58, 13), {date:'2020-10-04', count:1}),
  ee.Feature(ee.Geometry.Point(77.58, 13), {date:'2020-10-02', count:3}),
  ee.Feature(ee.Geometry.Point(77.58, 13), {date:'2020-10-01', count:2})
]

Should become this:
[
  ee.Feature(ee.Geometry.Point(77.58, 13), {date:'2020-10-02', count:4}),
  ee.Feature(ee.Geometry.Point(77.58, 13), {date:'2020-10-04', count:1}),
  ee.Feature(ee.Geometry.Point(77.58, 13), {date:'2020-10-01', count:2})
]

I have tried everything on planet Earth. I tried using dictionaries also, but when I try to check if a date is present in two dictionaries or not using .contains() and && operator it doesn't produce Boolean values correctly.


texts: PostGIS: Geography to geometry, best projection and architecture

texts:
Problem

I'm developing an app, that needs to calculate distances between my current position and different places in earth. I use PostGIS.

As my first solution I've stored all the positions as Geography, and I calculate distances easily. Everything is correct and working.
However, I need to do this queries lots of times, so I want to speed up this calculation.

Premises

I have two premises:

I store places around the world
I only need to check distances between my current location and places that are near (100km  distance max).
some other variables

We may have about 100-1000 places that match the distance
places do not move

Proposed solution

My idea is to use geometry instead of Geography.

I was thinking on using UTM.
For each place I'll store its position stored in its UTM zone projection

Each place that is 100km near the border of the zone, I also calculate its position in the neighbor zone and store it (This may put some error in the distance, but I don't think will be much).

When I want to look for distances I'll only check for the points that have the same zone, and calculate the distance.

As an example, I'll have a table like

| id    | place_id  | point     | utm_zone  |
|----   |---------- |-------    |---------- |
| 1     | 1         | 1,1       | 24M       |
| 2     | 2         | 0,0       | 25M       |
| 3     | 2         | -1,-1     | 24M       |

Question

Is this idea correct or I'm missing something?
Is there any other better way to do this?


texts: Convert shapefile to multiple single-value TIFFs

texts:
I have a shapefile that is composed of multiple overlapping polygons that represent extent of habitats of certain species. Here is a 1 MB example of such file from International Union for the Conservation of Nature website. I would like to convert the shapefile into multiple GeoTIFFs such that each file contains range map of a single species. How can I do that?

In ArcGIS Info I tried "polygon to raster", but the resulting raster is composed of a single band and multiple values (so that pixel values range from 0 to N, where N is number of species). This way information about overlap is lost. If I could convert my shapefile to multiband (such that each band represents a single species), the perhaps I could use the solution from this post. It would be best if my output images had 1 bit per pixel, since the range map should be essentially composed of zeros (species not present) and ones (species present).

Unless you know better options than ArcGIS? I also have at my disposal QGIS, GRASS and command-line GDAL tools.


texts: Shapefile empty after export of Add XY Data using ArcMap

texts:
I am using ArcMap.
I have a sheet with 3111 records. Using "add xy data" its attribute table is fine but when I exported it as .shp file it only displays header.
How do I solve this problem?


texts: How to query latitudinal bounds on 95% of a polygon in PostGIS?

texts:
I'm going to say thanks in advance. I'm very new to GIS and this forum is a tremendous help. My simple questions have been solved; here's one that might be a little trickier.

Given a polygon in a PostGIS geometry table, I want to find the northern and southern latitudes that would bound an arbitrary portion of the polygon's area (say 95%), with the remainder trimmed equally from the north and the south. I think what I want may be called a bounding box, though for a large polygon it might seem more like a conic section. I've noticed a number of functions relating to boxes in the documentation, but haven't figured out how they could be defined as I've described.


texts: Joining nearest points to make line in QGIS?

texts:
I have a dataset containing nearly 200 points.

Is there a tool for joining nearest points to make a multiline in QGIS?


texts: Lat/Lon point transformation/offset on arbitrary segment

texts:
I'm looking for something very similar to those stack exchange GIS questions:
Algorithm for offsetting a latitude/longitude by some amount of meters, How to convert (geodetic pose + cartesian offset) to geodetic position, Increment GPS coordinates by a certain distance given in meters [duplicate] and Quick way to determine if facing a given lat/lon pair with a heading, but I cannot get my head around to combine them for my specific question:

The most relevant would be the first one, but it refers to a specific vector/heading offset. What I'm looking for is a way to calculate an offset of some meters (less than 200 meters) for a point A, on the segment defined by this point A and another point B, having less than 1km of distance between them.

How can I apply the first answer of SE GIS with a calculation of heading (like the one in the last answer) in order to transpose point A for 100 meters in the direction of B?


texts: Computing distance raster in Google Earth Engine

texts:
I am trying to compute a distance raster in GEE. I am looking for a tool similar to Arcmap's 'euclidian distance'
I tried to use ee.Image.distance, however, I get an error related to the kernel.
Here is the the code: https://code.earthengine.google.com/5211d228ce97e647f9078eee7ad45ee2
This is the type of raster that I am looking for:

Thank you!


texts:
I'm trying to fix the order of appereance of my GeoJSON layers. Precisely, I have a layer of circle markers that I want always to be over polygon layers. At the moment, when a polygon layer is checked in the controller after the circle markers layer, it appears above... Does someone have an idea ?
Here is my layers control :
L.control.layers(baselayers, data, {position: 'topright', collapsed: true}).addTo(map);


texts: Identify features with temporal controller

texts:
I'm loving the new Temporal Controller!
One issue I have found is that the 'Identify Features' function does not seem to respect the temporal range specified in the temporal controller. In contrast, if I use a Filter (rather than the temporal controller) to restrict the displayed features then 'Identify Features' only selects displayed features.
Worth writing a bug report?
Edit in response to the comment requesting more detail follows:
I'm using QGIS 3.14.1-Pi from the QGIS Debian repository under Debian 10.
Details of the problem:
I have a temporal (single date column) PostGIS layer consisting of Polygons and MultiPolygons. I select a fixed range for temporal navigation and narrow my view down to only two polygons. I use the Identify Feature function to select one of those two polygons. QGIS responds with five different polygons, the one I wanted plus four more than are not within the fixed temporal range.
If it's really useful I can work on producing a minimal test case to demonstrate the issue. For now I'm hoping/assuming this is sufficient info to identify the problem.


texts: How do I put infowindow contents into a separate `div`?

texts:
I'm new to CartoJS. When a user clicks a dot on my CartoDB map, I want the content of the infowindow to appear in a separate div on the web page.

How do I do this with CartoDB/JS?